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ABSTRACT 



Traditional methods of time difference of arrival (TDOA) determination suffer 
significantly in environments fraught with co-channel interference and low signal to noise 
ratios. Cyclostationary signal processing techniques offer solutions to the shortcomings 
of the traditional TDOA methods. Specifically, the Spectral Coherence Alignment 
(SPECCOA) method of TDOA determination, developed by the Mission Research Corp. 
and Statistical Signal Processing Inc., performs exceptionally in very poor signal to noise 
ratio environments. The Applied Research Lab at the University of Texas at Austin 
(ARL:UT) has developed a prototype TDOA system, the Carry-on Multi-platform GPS 
Assisted Time Difference of Arrival System for the Naval Information Warfare Activity. 

It currently utilizes a traditional complex ambiguity function (CAF) to determine the 
TDOA(s) between multiple observers and an ARL:UT developed closed form solution for 
the geolocation of the emitter. The work presented here takes the first step in applying 
SPECCOA to the ARL:UT system. Coding both SPECCOA and the ARL:UT closed 
form solution in MATLAB® makes possible a quantitative comparison between the CAF 
and SPECCOA using ARL:UT real world test signals. 
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I. 



INTRODUCTION 



A. BACKGROUND 

The passive geolocation of electromagnetic (EM) emitters plays an increasingly 
important role in all three aspects of Information Warfare (IW). In order to protect 
friendly communications from hostile jamming and interference, information warriors 
must first locate the source of that jamming and interference. The geolocation of an 
enemy radar exploits the enemy by gathering information ascertained from his own 
information gathering system. Finally, in order to attack an enemy communications node, 
its coordinates must be passed to the strike aircraft or the cruise missile. In addition to its 
importance in the military arena of IW, passive geolocation of EM emitters finds use in 
law enforcement surveillance, search and rescue operations, navigation and the 
enforcement of communications regulations. 

Three aspects of EM waves lend themselves to exploitation for geolocation 
purposes. The path of an EM wave can be closely approximately from transmitter to 
receiver given the frequency of the transmission and some basic characteristics about the 
atmosphere through which it traveled. In addition, when EM waves arrive at two moving, 
spatially separated receivers, the receivers measure different Doppler shifts in the 
frequency of the transmission. Finally, the time that an EM wave arrived at two spatially 
separate receivers yields valuable information. These three characteristics naturally 
enough, are the basis for the three basic methods of geolocation, Angle of Arrival (AOA), 
Frequency Difference of Arrival (FDOA) and Time Difference of Arrival (TDOA). 

AOA geolocation involves the use of highly directional antenna arrays. By 
measuring the phase difference in the EM wave at each antenna element the receiver can 
calculate the direction from which the wave emanated. Drawing a line from the precise 
position of the receiver in this direction yields a vector containing the position of the 
transmitter called a line of bearing (LOB). Moving the receiver and taking another 
measurement or using an additional receiver located elsewhere produces another LOB. 

At the intersection the these LOB(s) is the transmitter. AOA is the oldest of the three 
methods and has been known by many names throughout the years of its use. Beginning 



just prior to World War II, the name most commonly associated with AOA was 
triangulation. With small measurement errors, the three vectors had a finite width when 
drawn and thus intersected to form a small triangle within which the emitter was located. 
This point is illustrated in Figure 1-1 below. 




Figure 1-1 - Angle of Arrival geolocation 

For moving collectors, the Doppler effect ensures that the receiver will measure a 
different frequency of arrival than it would as a stationary collector. Given two collectors 
that measure two different frequencies of arrival (FOA), the difference of these FOAs 
(FDOA) can be used to determine the position of the emitter. FDOA produces a locus of 
points called an isodop along which the emitter lies. This isodop represents all the 
locations from which the EM wave could emanate and produce the difference in Doppler 
shifts measured between the two receivers. Figure 1-2 shows a two dimensional 
depiction of a family of FDOA contours in the plane of the paper. Each contour 
represents a specific FDOA between the observers. The emitter could lie anywhere along 
the appropriate contour, and thus multiple measurements are required to determine its 
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location precisely. While these contours can be approximated as two dimensional when 
the emitter is on the surface of the earth, for an emitter in three dimensions they would be 




Figure 1-2 - FDOA contours 

surfaces. A receiver wishing to use Doppler for geolocation purposes must possess the 
ability to measure frequency more precisely than the smallest Doppler shift expected. 
Otherwise, the small Doppler shifts go undetected as they are below the noise in the 
frequency measurement. This has been the limiting factor in the usage of FDOA for 
precise geolocation. Measurement of frequency precisely is much more difficult and 
costly than the measurement of time precisely. 

The precise measurement of time has been revolutionized by the high speed 
digital technology of the Global Positioning System (GPS) and synchronization circuits 
for atomic time standards. Using the difference in arrival times of an EM wave at two 
spatially separated receivers, a locus of points called an isochron includes all possible 
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locations from which the EM wave could emanate and arrive at the two receivers at the 
two different times measured. In two dimensions, like the situation that approximately 



s(t) + n,(t) 

Observer #1 ► 



s(t - D) + n 2 (t) 
Observer #2 > 



TDOA Estimator 



Estimate 
of D 



Figure 1 -3 - Generic TDOA estimator 

exists when an emitter is constrained to the surface of the earth, the isochron can 
be represented by a hyperbola. In three dimensions, the hyperbola becomes a hyperboloid 
of revolution. A simple block diagram for a generic TDOA estimator appears in Figure 
1-3. 

The traditional methods used for TDOA determination include the general cross 
correlation (GCC) (in the absence of Doppler shift) and the cross ambiguity function 
(CAF) (both TDOA and FDOA simultaneously). These two methods perform well in 
benign signal environments. That is, with a signal of interest (SOI) with a high signal to 
noise ratio (SNR), in the absence of jammers and other signals not of interest (SNOI) in 
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the frequency band around the signal of interest (SOI), the GCC and CAF methods are 
superb estimators of TDOA. However, in a peacetime EM spectrum increasingly packed 
with signals, and a battlespace spectrum potentially flooded with SNOI(s), hostile 
jamming and weak SOI(s) the GCC and CAF methods suffer from severely decreased 
signal to noise and interference ratio (SNIR). In light of these weaknesses in hostile 
environments, a relatively new concept in signal processing called cyclostationarity has 
been applied to the TDOA problem. 

Cyclostationary processes are characterized by time-periodic statistics. A wide- 
sense cyclostationary signal specifically has a mean and an autocorrelation that vary 
periodically in time. These periodicities are usually hidden in traditional second order 
stationary processing techniques such as the autocorrelation function or its frequency 
domain counterpart, the spectral density function. However, when second order 
cyclostationary techniques such as the cyclic autocorrelation function and the cyclic 
spectral density function are used, wide sense cyclostationary signals readily display these 
hidden periodicities. These periodicities now revealed may be utilized for the detection, 
classification and parameter estimation of the signal. Among those parameters that may 
be estimated is the TDOA of a signal impinging upon two spatially separated antennas. 

The robustness of the cyclostationary technique for TDOA lies within the fact that 
every signal has a unique set of periodicities that depend on such parameters as 
modulation scheme, bit rates and spreading codes. Two signals that may appear as 
spectrally overlapping in the traditional stationary spectral density function can be 
separated with great accuracy by the cyclic autocorrelation function. In highly corrupt 
environments, cyclostationary processing techniques provide the capability to select 
specific signals for geolocation, nearly independent of the presence of jamming, inference 
and SNOI(s) that may spectrally overlap the SOI. 

The United States Navy has historically relied upon AOA techniques for direction 
finding. Beginning in the 1950’s, the Navy began constructing the High Frequency 
Direction Finding (HFDF) network of circularly disposed antenna arrays (CDAA). These 
CDAA(s) consist of a circular array of elements and reflectors that serve to detect the 
direction of incoming energy. Given the mass of recent base closures, their numbers have 
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dwindled significantly. In addition, they are not an organic asset for battle group 
commanders considering the current tactical communications environment that relies less 
on long-haul HF and more on Very High Frequency (VHF) and Ultra High 
Frequency(UHF) point-to-point links and Super High Frequency (SHF) satellite 
communications links. 

In addition to the HFDF network, the Navy added the OUTBOARD and 
OUTBOARD II systems to some of its Spruance-class destroyers. These too are 
primarily HFDF assets and have been followed by the current construction of COMBAT 
DF equipped Essex-class amphibious ships and Arleigh Burke-class guided missile 
destroyers. These three shipboard systems have proven extremely useful during the past 
years but as they have become the object of intense scrutiny with shrinking ship-building 
budgets and increasing decommissionings. Consequently, the Naval Security Group 
Support Activity (NSGSA), in June 1993, contracted with the Applied Research Lab at 
the University of Texas at Austin (ARL.UT) to develop an affordable, low-risk 
TDOA/FDOA geolocation system using commercial off the shelf (COTS) and 
government off the shelf (GOTS) technologies. 

The Carry-on Multi-platform Global Positioning System (GPS) Assisted TDOA 
System was tested for the first time in August 1995 after fifteen months of research and 
development. In further testing off the coast of San Diego, California, the prototype 
system geolocated HF, VHF and UHF emitters with a mean square error of approximately 
100 meters, using a CAF based TDOA determination algorithm. 

As noted above, the CAF has performed exceptionally well in the relatively 
favorable conditions during prototype testing but in theory, is susceptible to low SNR 
conditions and co-channel interference. Cyclostationary processing techniques for TDOA 
determination offer potential improvement to the performance of the ARL.UT system in 
highly corrupt environments. The evaluation of a cyclostationary TDOA algorithm in 
conjunction with the closed form geolocation algorithm used in the ARL:UT prototype is 
the primary objective of the work presented here. 

As a secondary objective, the cyclostationary TDOA determination algorithm and 
the closed form geolocation algorithm will also serve as the first two blocks for testing in 
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a developing geolocation software workbench at Naval Postgraduate School. The intent 
is to model the system in the Simulink® environment developed by The Mathworks Inc. 
and allow different algorithms to be exchanged by merely interchanging the appropriate 
Simulink® blocks. This will allow for the testing of all aspects of the geolocation 
problem from data conversion, filtering, AOA, FDOA and TDOA determination and 
various geolocation solutions. 

B. ORGANIZATION 

The traditional TDOA determination algorithms, the GCC and the CAF along 
with their disadvantages in corrupt environments are present in Chapter II. Following the 
illustration of the vulnerabilities of these two algorithms, the theory of cyclostationarity 
■ and TDOA determination by cyclostationary techniques are presented in Chapter III. 
Chapter IV depicts the closed form solution to the TDOA geolocation problem developed 
by Dr. Petre Rusu at ARL:UT for the prototype geolocation system which is described in 
Chapter V. 

Chapters VI, VII, VIII and IX constitute the bulk of the investigation, 
summarizing the new algorithms developed in MATLAB®, the test plan and the results 
of those tests including the plans for Simulink® block development. Finally, Chapter X 
draws some conclusions, offers explanations for the anomalies encountered and suggests 
areas for further research. 
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II. TRADITIONAL TDOA DETERMINATION 



A. GENERALIZED MODEL 

For the development of the two traditional methods of TDOA determination, 
assume a stationary signal from a remote emitter impinges upon two spatially separated 
antenna elements as illustrated in Figure 2-1. In the presence of additive white Gaussian 
noise (AWGN) noise, the two signals at receivers 1 and 2 can be modeled as 



where ni(t) and n 2 (t) are assumed to contain only AWGN with no strong SNOI(s) in the 
frequency band of interest for the SOI, s(t). A represents the complex valued relative 
magnitude and phase mismatch between the receivers; D is the TDOA of the SOI 
between the signals. 

It follows that the autocorrelation and cross-correlation functions are given by 



x(t) = s(t) + n,(t) 



( 2 - 1 ) 



and 



y(t) = A s(t- D) + n 2 (t) 



( 2 - 2 ) 



R x (z) = R s (T) + R nt (r) 



(2-3) 



^,(T) = |A| 2 ^(T) + i?„ 2 (T) 



(2-4) 



R yx (x) = AR x (x-D) + R ninz (x) 



(2-5) 



and the respective spectral density functions are 

S x (f) = S s (f) + S n> (f) 



( 2 - 6 ) 



S y (f) = \A\ 2 S x (f) + S ni (f) 



(2-7) 
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+ S n t n 2 (/ ) 



( 2 - 8 ) 



S JX (f) = A-S,(f)-e 



-72 n/D 



These relationships contain the parameter of interest, specifically, D, the TDOA. 
The next two sections show two different methods used to extract D from these 
equations. The first, the General Cross Correlation (GCC) function, can be used in the 
absence of relative motion between the two receivers and the emitter. The second, the 
Complex Ambiguity Function (CAF), can be used to estimate FDOA and TDOA 
simultaneously as is necessary when relative motion between the receivers and emitter 
exists. 




Of key importance to the development illustrated here is the statistical 
independence of s(t), n\{t) and ri2(t) and the absence of in band interference. While it may 
be argued that the signals and noise measured at the two receivers can be correlated to 
some extent, that scenario greatly complicates this case and is beyond the scope of this 



10 



demonstration. The derivations that follow are drawn from [1], [2], which treat the GCC 
more specifically and [3] and [4] which address the CAF more so than the GCC. All 
assume the three components to be uncorrelated. 

B. GENERALIZED CROSS CORRELATION 

Applying (2-5) in the absence of SNOI in the SOI band of interest and statistically 
independent s(t), nj{t ) and n 2 (t) yields 

R yx (T)=AR s (z-D) + R nin2 (z) (2-9) 

This function will peak at X- D, the TDOA between receivers 1 and 2. Because nj(t) and 
n 2 (t) contain merely AWGN, their cross-correlation term in (2-9) above reduces the SNR 
of the measurement but does not add interference in the form of SNOI(s). Further 
analysis using (2-6) - (2-8) in a similar manner reveals 



S x (f) = S s (f) + S ni (f) 


(2-10) 


S > (f) = \A\ 2 S s (f) + S n2 (f) 


(2-11) 


S yx (f) = A-S s (f)-e- i2 ^+S nin2 (f) 


(2-12) 


To extract D from (2-12), first note 




S s (f) = < >0 

[=0 otherwise 


(2-13) 



which assumes the SOI effectively occupies a finite bandwidth B around the carrier or 
intermediate frequency / 0 . Taking the ratio of the spectra and using (2-10) with (2-13), 
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(2-14) 



y/) 

s x (f) 



A e j2vfD 

0 



/ - — </</ + — 
^ o 2 J J » 2 

otherwise 



and taking the inverse Fourier transform of this ratio gives 



f + B A 

J ° /2 C ( f\ 

/i„(x)= | y- e +i2Kf, df 



f>- B A 



s x (f ) 



A sin[7i B(x- D)] 
K(x- D)/2 



cos[27t/„(T-D)] 



(2-15) 



which peaks at T= D. This in turn may be rewritten as 



l + b A 

KA) = JW)^(/)^# (2-16) 



The weighting function, W(/), is defined in this instance as 1 /£*(/) in (2-16) which 
is the best choice given no a priori information about the SOI [5]. In addition, it 
distinguishes this case as the generalized cross correlation method. Assigning a 
weighting function, W(J) = 1 reduces the TDOA determination to a simple cross- 
correlation as in (2-5). Given prior knowledge of the noise and interference 
characteristics, other choices for W(f) include the Roth impulse, SCOT and PHAT among 
many others [3] which are designed to reduce specific noise and interference problems. 

It is clear that, in the absence of significant noise, co-channel interference and 
relative motion between receivers and emitter, the GCC produces the desired estimate of 
the TDOA of a signal from the ratio of the estimates of the spectral density function of 
the signal at one receiver and the cross spectral density function of the two received 
signals. Figure 2-2 above illustrates the generalized cross correlation process. The two 
filters, H\(f) and H 2 (f) are specifically designed to remove out of band interference. 
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While in practice the GCC is often done in the frequency domain. Figure 2-2 depicts the 
time domain correlation of the input signals. The processes are theoretically equivalent. 

Clearly, the presence of an interference signal in the bandwidth B defining the 
spectral densities would corrupt the estimate. In addition, as is shown in the next section, 
Doppler shifts at each receiver also must be accounted for in order to accurately 
determine TDOA. 




Estimate of TDOA 



Figure 2-2 - GCC block diagram 



C. COMPLEX AMBIGUITY FUNCTION 

The complex ambiguity function (CAF) can be interpreted as an extension of the 
GCC for moving transmitters and/or receivers. Stein in [4], showed that in order to 
properly determine the TDOA between two receivers in the presence of a Doppler 
difference at each receiver, the spectrum of one of the signals first must be translated in 
frequency by an amount/equal to the Doppler difference (FDOA) measured between the 
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observers. In order to show this, a Doppler shift is introduced into the generalized model 
from above. Equations (2-1) and (2-2) can now be written as 



40 = s (0 + > i ,(0 



(2-17) 



y(t) = As(t - D)e~ J2l9 ‘ 1 ' +n 2 (t ) 



(2-18) 



where fd is the Doppler difference as measured between observer 1 and 2. 

Now, in place of calculating TDOA with the two dimensional GCC from (2-16), 
which in the presence of significant Doppler could peak at a value that does not truly 
correspond to the TDOA, it is necessary to calculate the three dimensional CAF given by 



The simultaneous determination of the TDOA D and the Doppler shift /d causes 
\A(D,f d )\ to peak. At/d = 0, the CAF reduces to a GCC problem as outlined above. For 

fd* 0, the CAF may be thought of as a GCC performed after frequency shifting the 
spectrum of y(t) up or down as necessary by an amount equal to //. A block diagram of 
the CAF operation appears in Figure 2-3. Note the similarity between it and the GCC 
diagram of Figure 2- 1 . 

The three-dimensional width of the correlation lobe is directly proportional to the 
accuracy of the estimates of TDOA and FDOA. Stein further points out in [4] that the 
variance of the estimate for each parameter can be related to the noise bandwidth, the 
signal bandwidth, the integration time and the effective input signal to noise ratio as 
follows 



T 




(2-19) 



0 



PV57/ 



( 2 - 20 ) 
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Figure 2-3 - CAF block diagram 



where 



1 1 

aFDoA " x 7m 

B = noise bandwidth at the input of both receivers 



(5 = 2 it 



\f 2 W s (f)df 

\W s (f)df 



Vi 



where W S (J) is the spectral density of the signal as shaped by the receiver. 

T e = rms integration time 



( 2 - 21 ) 



( 2 - 22 ) 
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(2-23) 



--- 1 , 1 | 1 
Y 2 |_y, y 2 y,Y 2 _ 

where — and — are the SNR(s) for each receiver respectively. 

Y. y 2 

D. PERFORMANCE 

Given (2-20) and (2-21), it is clear that both the GCC and the CAF can be 
rendered ineffective by moderate in-band noise. Figure 2-4 below is an illustration of the 
theoretical standard deviation of a TDOA measurement made at various SNR(s) and 
integration times. It considers only the errors introduced by AWGN according to the 
theory presented above and includes no error from any other source such as 
instrumentation or machine precision for example. Because theory dictates these values 
as the minimum errors, they may be considered the lower bound for this scenario. 

The model for this demonstration assumes a rectangular signal spectrum with 
rectangular, full bandwidth Gaussian noise with no Doppler between observers. It 
assumes no magnitude or phase mismatch between receivers and the same SNR at both 
receivers. It approaches the ideal situation given the perfect match between signal 
bandwidth, noise bandwidth and receiver measurements. In reality, the signal bandwidth 
would be more Gaussian or raised cosine pulse shaped and the receivers’ measurements 
would not be matched in phase or magnitude. These differences would corrupt the 
TDOA estimation further, requiring more SNR and integration time to achieve the lower 
bounds determined by this idealized model. 

Clearly, below approximately 10 dB SNR and 400 ms the idealized model 
exceeds a TDOA standard deviation of 100 m. Given three TDOA measurements at this 
level and a simple Euclidean norm combination of these errors, a geolocation derived in 
the absence of any geometric dilution (see Chapter V) would have a standard deviation on 
the order of 173 m. This exceeds the 100 m goal of the ARL:UT project without 
considering the additional errors introduced in an actual system (again see Chapter V). 
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Figure 2-4 - Theoretical standard deviation of TDOA measurement given various 
SNR(s) and integration times 

Poor SNR(s) and SNIR(s) are obviously problematic for the GCC and CAF given 
the need for geolocation accuracy on the order of 100 m. While increased integration 
times can solve some low SNR problems, this solution has an upper bound. The 
integration time must not exceed the coherence of the SOI. That is, the statistics of the 
SOI must remain relatively stationary during the integration time for the GCC and CAF to 
operate properly. A significant increase in integration time to combat poor SNR may 
exceed the coherence time of the SOI and introduce yet more error. Modeling the SOI as 
cyclostationary vice stationary and employing the appropriate processing techniques can 
in many cases overcome most of these problems. 

Cyclostationary techniques exploit periodicities introduced to man-made signals 
in a number of ways. These periodicities can unique to specific signals and thus can be 
used to distinguish one signal from another in the same bandwidth as well as significantly 
reduce the level of post-processing noise. Thus, with cyclostationary signal processing, it 
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is possible to tolerate significantly lower SNR(s) and still obtain excellent TDOA 
measurements. Chapter III introduces the theory of cyclostationarity in communications 
signals and develops a method for TDOA determination. 
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III. CYCLOSTATIONARY TDOA DETERMINATION 



A. DEFINITIONS 

The concepts of cyclostationarity have been examined in theory for over two 
decades. Beginning in the late 1960’s, Dr. Franks of the University of Massachusetts at 
Amhurst and Dr. Gardner of the University of California at Davis (UCD) began extensive 
research in the area of cyclostationary signal processing. Dr. Gardner has since produced 
multiple publications in the field. His paper on the general theory of cyclostationary 
signal processing, published in the April 1991 edition of the IEEE Signal Processing 
Magazine [6], stands as the key reference for most aspects of the theory. As is always the 
case when one individual has contributed so much to an important engineering 
development, much of the theory and examples presented here follow Dr. Gardner’s work 
closely. His solo efforts and collaboration with others appear in references [6], [7], [8] 
and [9] and is redundant in many instances. Where appropriate, the specific source of 
unique information is referenced below. 

As previously noted, most modern signal processing techniques associated with 
communications applications treat SOI(s) as stationary random processes. However, 
because most manmade signals are generated through some repetitive, periodic process 
such as the amplitude, frequency or phase modulation of a sinusoidal carrier, the 
encoding of data or the encryption of a message, their statistics inevitably vary 
periodically with time. While in many instances, receivers may successfully ignore the 
underlying periodicity of a manmade signal, often the detection of the signal and the 
estimation of its parameters is more successfully accomplished by modeling the signal as 
cyclostationary vice stationary. 

Simply stated, a process with statistics that vary periodically with time is termed 
cyclostationary. Figure 3-1 depicts a block diagram of the procedure that leads to a 
cyclostationary signal for most communications processes. A stationary random message 
such as digital data or analog voice is modulated, clocked or framed by a periodic 
communications process through some nonlinear system and a cyclostationary signal 
results. 
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Stationary 
Random Message 

- digital data 

- thermal noise 

- analog voice 

Cyclostationary 
Signal 

Source of Periodicity 

- RF oscillator 

- bit-stream clock 

- repeating frame structure 

- rotating machinery waveforms 

(motors, engines, turbines, propellers...) 




Figure 3-1 - Origins of cyclostationarity, adapted from [10] 

From a strictly mathematical point of view, a cyclostationary signal of order n is 
one that will have additive sine-wave components that result in spectral lines for some n Ih 
order nonlinear transformation of the signal. In the case of n = 2, a signal is said to be 
second order cyclostationary if a quadratic transformation produces additive sine-wave 
components that generate spectral lines. This characteristic may be thought of as akin to 
a process being considered wide sense stationary or stationary through order 2. 

To continue the case of n = 2 more specifically, a signal x (t) is cyclostationary 
with cycle frequency a if and only if some of its delay products y(t) = x(t)x(t - x) produce 
a spectral line at frequency a. If not all cycle frequencies a for which x(t ) exhibits 
cyclostationarity are harmonics of a single fundamental frequency, then jc(r) is 
polycyclostationary. Polycyclostationarity implies the existence of more than one 
periodicity in the statistics of a signal. This in turn implies more than one source of 
periodicity such as the case when clocked-digital-data phase modulates a sinusoidal 
carrier in an M-ary phase shift keying scheme [6]. 
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To further illustrate the concept, consider a signal x(t) that contains an additive 
sine-wave component at frequency a and is of the form 

a cos(2otf + 0) with a ^0 (3-1) 

The Fourier coefficient defined as 



M“ = (x(t)e~ i2KaJ ^ 



(3-2) 



where (•) denotes time-averaging, will be non-zero. Note that this is of the form of the 
common representation of the Power Spectral Density (PSD) of x(t ) with a spectral line at 
+a and its mirror at -a. In simple terms, x(t) contains first-order periodicity at frequency 
a given by 



| M°\ 2 [8(f-a)+5(f + a)] (3-3) 

Now, considering contributions to x(t) from other sources, the total signal may be 
represented as 



x(t) = a cos(27tar + 0) +n(t) (3-4) 

where n(t ) can be thought of as the random energy outside that of the SOI. If n(t ) is 
strong in comparison to the SOI such that it masks the sine-wave components in x{t) from 
detection during casual inspection of the waveform, then x(t) can be thought of as 
possessing hidden periodicity. This hidden periodicity can still be exploited by using 
signal processing techniques such as the PSD function as noted above. However more 
powerful cyclostationary signal processing techniques are more sophisticated and can 
unmask periodicities in signals that are hidden even from common traditional techniques 
like the PSD function. 
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B. CYCLIC AUTOCORRELATION FUNCTION 



Traditionally, a common second order, time-domain statistic used in signal 
processing is the autocorrelation function given by the quadratic transformation 

R xx (t:) = (x(t)x(t - r)) (3-5) 

In the case of a signal that can be modeled as cyclostationary, this transform will produce 
spectral lines at non-zero frequencies a such that 

M“ =(y(0e‘' 2 ™)* 0 (3-6) 

where 

y(t) = x(t)x(t-x) (3-7) 

This signal x(t) as noted above contains second order periodicity manifested in the 
PSD of the delay product given in equation (3-5) above. Transforming (3-5) into a 
symmetric delay product and accommodating complex signals as well gives 



y(0 = * 



f t) 


* 


f A 




X 




l 2 J 




l 2) 



(3-8) 



which forms the basis for the fundamental second order cyclic moment known as the 
cyclic autocorrelation function: 



KW 





(3-9) 



Of note, the cyclic autocorrelation function assumes the form of the Fourier 
coefficients of the additive sine-wave components produced by the periodicity of the 
delay product in (3-6). 
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Another interpretation of (3-9) is as the traditional stationary autocorrelation 
function multiplied by a kernel e i2KCa that produces spectral lines at frequencies where the 
stationary autocorrelation function contains additive sine-wave components indicative of 
its periodicity. Consequently, for either interpretation, at a = 0, the cyclic autocorrelation 
function reduces to the traditional autocorrelation function. 

A final interpretation of the cyclic autocorrelation function can be seen by 
defining first 



u{t) = x(t)e 



(3-10) 



and 



v(r) = x(t)e +]KOJ 



(3-11) 



so that u(t ) and v(?) are frequency shifted versions of x(t). Now /?“ (x) can be written 



K(r) = (u 




(3-12) 



which is the cross-correlation of the two versions of x(r) shifted up and down by 
frequency a/2. In other words, the cyclic autocorrelation function may be viewed as the 
correlation in the time domain between two values of x(t) separated in frequency by a. 
Consider as an example, a Binary Phase Shift Keyed (BPSK) signal given by 

s(t) = A c d(t)cos(2nf l t + <^>) (3-13) 



where d(t) is the binary modulating wave form consisting of positive and negative 
rectangular pulses given by 
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(3-14) 



d(t)= ^d„q(t-t 0 -nT b ) 
« = -«> 



Approximated as a random binary wave, it contains no first order periodicity and thus no 
spectral lines in its PSD. Consequently, the PSD of the BPSK signal 5(0 is a scaled sine 
squared function and will also contain no spectral lines as is evident its expression 



Sbpsk(/) = 4 2 7; 



r sin jtf c T b Y 









(3-15) 



Thus, the autocorrelation function of the BPSK signal, the inverse Fourier transform of 
(3-15), will contain no additive sinusoidal components to indicate any periodicity. 
Multiplying the conventional autocorrelation function by the cyclic kernel e'* 2nat , the 
cyclic autocorrelation function follows as 



K = ~ ''“('*") cos(2jtf c r)e~ j2mT ° 



L[ r «+2/c ( T ) e -> 2)r [“ +2 /c lT » _|_ r “-2 fc ^ e -j 2 ”l a ~ 2 fc^o j 



(3-16) 



47; 



where 







(3-17) 



This clearly shows additive sine-wave components at a = ±2/ f ± R h and a = ±R h [9]. 
Thus, the quadratic transform that is the cyclic autocorrelation function unmasks hidden 
periodicity in this simple BPSK signal and proves it polycyclostationary in the process. 
Naturally, the conventional cross-correlation function and the cyclic cross-correlation 
functions are related in a similar manner. 
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c. 



SPECTRAL-CORRELATION FUNCTION 



The frequency domain equivalent for the cyclic autocorrelation function is the 
spectral correlation function (SCF). It follows directly from the Fourier transform of (3- 
9) and has the form 



1 W+A//2 

Sx T (t,f) Al =-—\ X T (u,f+a/2)X T (u,f-a/2)du (3-18) 

where 



rt+T/2 ... 

X T (t,f)= I x(u)e- }2nfu du 

Jt-T/2 



(3-19) 



is the finite-time Fourier transform of x(t). More commonly used in cyclostationary 
signal processing, the SCF reduces to the conventional PSD at a = 0 just as the cyclic 
autocorrelation function reduced to the conventional autocorrelation function at a = 0. 

Continuing the BPSK example from above, the SCF can be found directly from 
(3-17) or by taking the Fourier transform of (3-15). It has the form 



S“(f) = 

* 471 



J r , 

Q f + fo+— Q /+/„-— 

V 2.) V 2) 



^ 0(,\ .( CC ^ 

+Q /-/„+— \Q f - f„~ — 

V 2) V 2) 



- j2Koa 0 



for a = n/T h and 



(3-20) 



W) = 



a 



4 T b 



+Q 



a 



Q { f + f ‘’ + 2 ) Q [ f ~ L ~2 



r , 0C\^( r r OC ^ 

f +fo + “|2 ~ f '~2 



-j2K[a+2f c ]t 0 



-jlK[a-2f c )t 0 



for a - ±2 f c ±n/T h with Q(f) being the Fourier transform of the keying envelope q(t ) 

[ 8 ]. 
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The SCF derives its name from the fact that it is in fact a correlation of the SOI in 
the frequency domain. The SCF at frequency / 0 and cyclic frequency do is merely the 
correlation of two values of the signal in the frequency domain separated in frequency by 
Oq and centered at frequency f 0 . Figure 3-2 shows this process concept graphically. 




Figure 3-2 - Spectral Correlation Function illustrated, from [10] 

The SCF is plotted on what is called the bi-frequency plane. The plane is defined 
along one axis as spectral frequency / and along the opposing axis as cyclic frequency a. 
Figure 3-3 illustrates this point. The magnitude of the SCF corresponds to a height above 
the bi-frequency plane and is often plotted in that fashion. Note that the values for a 
range from -f s to/ s while values for spectral frequency /naturally range from -// 2 to/ s /2. 
Because a corresponds to the separation distance of the correlated values in the frequency 
domain, its range extends twice that of/. 

The SCF computation can be highly complex and demanding on any processor. 
For that reason, two algorithms were developed in [7] and designed specifically for 
computation efficiency. The Fast Fourier Transform Accumulation Method (FAM) and 
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Figure 3-3 - Bi-frequency plane 
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Figure 3-4 - SCF of some common modulated signals, adapted from [10] 
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the Strip Spectral Correlation Analyzer (SSCA) both simplify and thereby accelerate the 
computation of the SCF. Because both the SCF and the Cross Spectral Correlation 
Function (CSCF) are necessary for cyclostationary TDOA algorithms, SSCA plays an 
important role in this work. It is presented in more detail in Chapter VI. 

Finally, in addition to their role in the TDOA computations, the SCF and CSCF 
can also be used for purposes ranging from signal detection to characterization and 
estimation of parameters. Specifically, the modulation type of an unknown signal can be 
found using the SCF of that unknown signal and comparing it to SCF(s) of known 
modulation type. Figure 3-4 shows a small library of SCF(s) from several more common 
modulation schemes. 

D. CYCLOSTATIONARY TDOA SIGNAL MODEL 

The signal model used to develop cyclostationary TDOA algorithms is very 
similar to that presented in Chapter II. Recall the signals received at two spatially 
separated observers can be given by (2-1) and (2-2) and are below for clarity 

x(t) = s(t)+n,(t) (3-21) 

and 

y(t) = A- s(t — D) + n 2 (t) ■ (3-22) 

The difference between the two models lies in the temporal and spectral relationships 
between s(t), tii(t) and n^t). In Chapter II, «i(r) and ni(i) were assumed to have no 
temporally and spectrally components that overlapped the SOI, s(t). In the 
cyclostationary TDOA signal model, n\(t) and «2(0 represent all SNOI(s) and noise 
present at the respective observers. They may or may not contain co-channel interferers 
in reality, however for purposes of this work, they are assumed to contain interference 
that in fact spectrally overlaps the SOI. Figure 3-5 portrays this situation graphically. 
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Figure 3-5 - Cyclostationary signal model scenario 



Given this model, the cyclic auto and cross-correlation functions in general can be 
written as 



/?;(T) = /?,“(T) + /?“(T) 


(3-23) 


^(T) = |A| 2 ^(ry- w + /?“(T) 


(3-24) 


(t) = ARf (T - D)e~ jKaD + R“ ni (r) 


(3-25) 



just as (2-3) - (2-5) represented the conventional auto and cross-correlation functions. In 
addition, as with (2-6) - (2-8), the spectral correlation functions are 

W) = S“(/) + S“(/) (3-26) 
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1 S“(/) = |A| 2 ^(/y 2 ”“ D +5“(/) 



(3-27) 



S; x (f) = AS:(f)e 



- j2jt(f +af2)D 




(3-28) 



Though (3-26) - (3-28) contain the TDOA of the signal D, this parameter is now 
masked by the spectrally overlapping components of ii\{t) and / 12 (f)- Thus, any attempt to 
estimate D with traditional methods that compute the above equations at a = 0, like those 
illustrated in Chapter II, would result in a corrupted value. However, if s(t ) contains 
some cycle frequency, oto, that it does not share with any component of n\(t) and 
then cyclostationary techniques can discriminate those contributions to D made by s(t ) 
from those contributions of n\{t) and / 12 (f) that would otherwise corrupt the estimate. 
Essentially this eliminates the SNOI(s). Thus a reliable estimate of D is possible even in 
highly corrupt environments provided a unique a exists for the SOI. Spectral Coherence 
Alignment is a cyclostationary TDOA algorithm employing measurements of both the 
SCF and CSCF. It determines the TDOA(s) for the TDOA processor developed here and 
appears in detail below. 

E. SPECTRAL COHERENCE ALIGNMENT 

Spectral Coherence Alignment (SPECCOA) was developed and presented in [8] 
using an ad hoc minimum least squares (MLS) approach. Using equations (3-23) - (3- 
25), and the assumption that s(t ) contains a unique a, it can be seen that 



The cross-correlation term for / 21 (f) and / 12 (f) is eliminated by the use of a cycle frequency 



unique only to the SOI. An estimate of the TDOA, D that minimizes the sum of the 
squares of error magnitudes between the measured value of the left side of (3-29), /?“ 



R a yx {u) = CR*(u-D)e- 1KaD 



(3-29) 
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and the measurement of the right side of (3-29), R“ with D substituted yields the MLS 
optimized value for the TDOA. 

Mathematically, from [8], 



A A 

D=argmin<{ c a ( t) 



(3-30) 



A 

where c«(r) is the estimate of the MLS function. It was further be shown in [8], that the 
solution of the MLS problem is of the form 



where 



A A f A 

D=argmaxj c a ( r) 



c’ a {T)=\\R a yx {u)R a x {u-T)'du 



(3-31) 



(3-32) 



= \js; x (f)s^(fYe j2 ^df\ 



(3-33) 



where (3-33) derives from (3-32) through Parseval’s relation. Gardner and Chen go on to 
further prove that (3-33) does indeed peak at r = D. 

SPECCOA derives its name from the fact that the peak in (3-33) occurs by 
maximizing the correlation in / of the two spectral correlation functions through the 
alignment of the phases of the respective functions. 

SPECCOA proves more useful in tactical applications than other cyclic TDOA 
algorithms (see [5] and [9]) by virtue of its ease of implementation and performance in 
corrupt environments [9]. Though other cyclostationary algorithms consistently out 
perform SPECCOA [5], they do so at the expense computational complexity and speed. 
Because this work is intended to demonstrate the feasibility of implementing a 
cyclostationary TDOA algorithm in an operational tactical system, SPECCOA was the 
logical choice for its efficiency and performance. Chapter IV discusses the algorithm 
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which utilizes TDOA(s) to determine the geolocation of an emitter. Chapter V presents 
the ARL:UT prototype TDOA geolocation system. Chapter VI contains specific 
implementation issues for SPECCOA in the MATLAB® environment. 



IV. TDOA GEOLOCATION CLOSED FORM SOLUTION 



A. BACKGROUND 

Given the TDOA(s) from pairs of receivers, the problem now lies in determining 
the location of the transmitter from the intersection of the isochrons generated. Many 
distinct methods appear from the geometric interpretation of Schau and Robinson in [1 1] 
to the iterative solution proposed by Loomis in [12] and the closed form solution 
illustrated by Ho and Chen in [13]. While somewhat similar, they all take a specific 
course in determining the location of the transmitter given the determined case. What 
becomes more difficult for each of these solutions is the determination of the contribution 
of errors in the measurement of the TDOA(s) to the final result, the geolocation. 

In [14], Rusu develops a closed form solution avoiding the initial guesswork 
involved in an iterative technique. Unlike [13], his solution is completely general. In 
addition, he shows that the errors in the measurements may be propagated through the 
mathematical model in a linear fashion, simplifying the determination of uncertainty in 
the geolocation. The following development draws exclusively from his derivation. 

The problem of determining the location of an emitter in three dimensions given 
the times of arrival (TOA) at four spatially separated receivers has often been treated as a 
TDOA problem. That is, time is treated as absolute time. However, as Rusu points outs, 
the mathematical model is invariant to time translation and thus may be treated as a TOA 
problem by setting the arrival time at any one of the receivers as the origin of the 
problem, t = 0. This simplification without loss of generality facilitates the development 
of linear error propagation from the initial measurements through the mathematical model 
to the final result as noted above. 

In the presence of moving emitters or receivers, both TDOA and FDOA must be 
determined by one process as shown in the development of the CAF in Chapter II. The 
case of FDOA can be similarly simplified by assigning a Doppler shift of zero to one of 
the receivers reducing the problem to one of frequencies of arrival (FOA). Fortunately, 
given the independence of the TOA portion of the model from the frequency of 
transmission and Doppler shifts, the TOA and FOA equations may be treated separately. 
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Though TDOA is the emphasis of this thesis, for completeness both the TOA and the 
FOA solutions will be developed here. 

B. MATHEMATICAL MODEL 

The explicit functional model relating the observable X to the dependent variable 
of interest Y, in the TOA case is the solution of the set of scalar equations of the form 
F(X,Y) = 0. F may be referred to as the implicit functional model, X the independent 
variable and Y the dependent variable. In the development of the TOA-FOA problem, X 
is a 32-dimensional real vector state space of observations while Y is an 8-dimensional 
status vector state space. Figure 4-1 illustrates the scenario in general. 




The TOA-FOA case assumes that the i-th receiver located at rj = (Xj, yj, Zj) and 
moving with velocity v \ = (Uj, Vj, wO detects a signal at time r, with frequency co,. Four 
observers for the determined 3-dimensional case, constitute the observation state space 
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vector of 32-dimensions. The parameters of the transmitter, its location and velocity 
given by r = (x, y, z) and v = («, v, vv) respectively, its transmit time t and frequency co 
compose the 8-dimesional status vector state space. 

The fundamental equation relating the signal travel time to the separation distance 
between observer i and an emitter assuming the signal travels at the speed of light is 

V (*.- “ • x ) 2 + O '.- “ y ) 2 + ( Z / _ -) 2 = c ( ? , “ 0 ( 4 - 1 ) 

where ■ y j(x i - x) 2 + (y,- - y) 2 + (z,. - z) 2 is the Euclidean distance between r n and r denoted 
hereafter by ||r n — rj| and c is the speed of light. Forming appropriate vectors using a 
single equation to include all four observers produces the implicit functional TOA model 

F; (X, Y) = ||r, - r|| - c{t, - 1 ), / = 1, 2, 3,4. (4-2) 

Simple algebraic transformation leads to 

\_ 

[(jc ; -A r) 2 + (y,. -y) 2 +(Zj -z) 2 -(t, -t) 2 ] 2 =0, / = 1,2,3, 4. (4-3) 



In vector form the TOA equations are defined as F T = (F\, Fj, F$, F 4 ). 

The Doppler equations can be similarly treated. The fundamental relationship 
between the Doppler shift observed by receiver n and the relative velocities of the emitter 
and observer is 







c J 




(4-4) 



where G>d is the Doppler shift observed, co is the transmitted frequency, (v n - v) is the 



r — r 

relative velocity and j- 2 r is the unit vector from the observer to the transmitter. 

r_ - r 
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Again forming vectors based upon this relationship yields the implicit FOA functional 
model 



F i+4 (X,Y) = co i -(0 



l-(v.-v) -5-1 
t. -t 



i = 1,2, 3, 4. 



(4-5) 



defined as F f = (F 5 , F 6 , F 7 , F 8 ). 

Combining the TOA and FOA models produces the TOA-FOA implicit functional 
model defined as F = (Ft, F f ). 

C. DETERMINATION OF THE TRANSMITTER STATE VECTOR 

The determination of the 8-dimensional transmitter state vector involves finding 
for X, the parameter Y that satisfies the combined TOA-FOA implicit functional model, 
F,(X,Y) = 0, i = 1,2. .., 8 , where the functions F, are defined in (4-3) and (4-5). The 
solution of such a system of equations occurs in two steps. First, the TOA equations can 
be solved for Y T = (r, t ). The resulting solution can then be used in the FOA equations to 
solve for the unknowns Y F = (v, co). 

The TOA equations are irrational in their form in (4-3) and thus must be squared 
in order to be solved. This produces a set of quadratics for which two sets of solutions 
exist corresponding to the two roots of the equations. In most cases, Rusu points out, the 
two solutions to the rational quadratics are also solutions to the original irrational set of 
equations. Rarely, most often in the 2-dimensional case, only one of the two solutions 
also solves the irrational set of equations and leads to a unique solution to the geolocation 
problem. 

Since each solution to the TOA case also produces a unique solution to the FOA 
equations, more often than not, there are two distinct solutions to the problem. This is the 
classic case of ambiguity encountered with every TDOA solution proposed thus far. 

Rusu handles this problem by noting that information outside the algorithm must resolve 
the ambiguity. Prior knowledge of the general area of the target’s position might 
eliminate one of the solutions. Multiple measurements may also reveal one solution 
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converging in one location while the other solution diverges and produces seemingly 
unrelated geolocations. 

To solve the TOA equations, (4-3) must be squared to produce 

(A - x) 2 + (y, - y) 2 + (z, - z) 2 - (t, - 1 ) 2 =0, i = 1, 2,3,4. (4-6) 

The solution to (4-6) is found by subtracting the equations two-by-two to eliminate x 2 , y 2 , 
z and r. This produces a set of linear equations in r and t that may be written as 

Ar = qr + s (4-7) 



As the equations can be subtracted in any order to produce such an elimination, 
there exists many possible resulting systems of equations. One more obvious possibility 
used by Rusu is 



A = 



x i — x 2 y, - y 2 z,-z 2 ^ 
X 2 x s y 2 -ys z 3 -z 3 
x 3 X 4 Vi - Vi z 3 -z 4J 



q = 



f t — t ^ 
l \ l 2 

t 2 — 1 3 

\h~Uj 



s = — 
2 



m 2 it ||2 2 

6 1 -hi ~ t i +t 



r 2 ~ H 



-nr, 



h + h 

— t 2 +t 2 
h — *4 



(4-8) 



(4-9) 



(4-10) 



This system makes it possible to solve for r in terms of t producing a t- 
parametrized solution to the TOA equations 

r = gt + h, (4-11) 
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where 



g = /I 1 q and h = A‘'s. 



(4-12) 



The introduction of this result into the k lh range equation from (4-6) yields a 



quadratic in time t. This equation, with some algebraic manipulation can be written as 



The two roots that are the solution to (4-13) produce two values for r when 
inserted into (4-1 1). Rusu notes that the roots are usually real but in the rare cases when 
they are complex, experience has shown him that the real part of the complex roots 
suffices as a solution. 

Once the TOA equations provide the emitter position and time of transmission, 
these values produce the transmitter’s velocity and transmitting frequency from the FOA 
equations. If the FOA equations are rewritten as 



at 2 + 2bt + c = 0 



(4-13) 



where. 




(4-14) 



b = g-h + g-r k +f t 



(4-15) 




(4-16) 




(4-17) 
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the FOA case becomes a non-homogeneous linear system that may be solved using 
standards procedures like Gauss-Jordan elimination [14]. 



D. ERROR PROPAGATION 

Perhaps the most valuable portion of Rusu’s solution comes with the error 
propagation. He developed a linear method through which errors in the measurement can 
be propagated through the model and result in predictable errors in the fix. The following 
derivation is based solely on his development of this method. Though he develops both 
the TOA and the FOA propagations, only the TOA case will be presented here. 

The key to the linear propagation of errors is based upon the assumptions. First, 

Y = Y(X) is assumed to be differentiable in the neighborhood of the measured value X 
containing the exact (error free) value of the observable X. Second, if the errors in the 
observation and the parameter are be denoted 5X = X - X and 5Y = Y - Y, Y is assumed 
to be linear in the region of these errors and thus 8Y can be written 



8Y = 



^dY 

v dX 



5X 



(4-18) 



The following Jacobian matrix represents the derivative of a vector valued 
function, Y = (Y|,}' 2 ,--.,l'n) with respect to a vector variable X = (;ti,X 2 ,...,*n) 



UxJ 



fd d 
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•w* J 



(4-19) 



Given the linear approximation in (4-18), the covariance matrix of the parameter 
can be related to the covariance matrix of the observation through 
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(4-20) 



r 

y laxj \axj 

The implicit function theorem provides a convenient method for computing the 
derivatives of Y with respect to X from the relation F(X,Y) = 0. This provides the 
general formula 



P Y 1 
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UxJ 


UyJ 


UxJ 



(4-21) 



Just as the relationship of the TOA and FOA solutions allowed the separation of 
the TOA case from the FOA case, the matrices involved in (4-21) can be written such that 
the TOA and FOA contributions are separated as 
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(4-22) 



(4-23) 



(4-24) 



The independence of the TOA solution from the FOA equations is evident in the 



zeros found above corresponding to 
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The problem of inverting the Jacobian, 3F/3Y, can be solved numerically 
however, Rusu managed to express the inverse of 3F/3Y in terms of the inverses of the 
TOA-FOA blocks forming it. He points out that this formally separates the TOA and 
FOA error propagations and allows computation of just the TOA solution and errors if 
just the transmitter location is desired. The block inverse of (4-24) is 
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(4-25) 



It is now possible to write ffie components of the right hand side of (4-22) as 
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(4-29) 



using (4-21), (4-23) and (4-25). 

Now the TOA Jacobian matrices related to the TOA error propagation can be 
explicitly expressed using the functional model given in (4-3). The Jacobians can be 
expressed as 
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where 



Implementation of Dr. Rusu’s algorithm into MATLAB® followed naturally 
given its matrix and vector nature. It plays an integral part in the ARL:UT prototype 
TDOA geolocations system which is discussed in Chapter V. Specifics regarding the 
encoding of Dr. Rusu’s algorithm can be found in Chapter VI. 
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V. THE CARRY-ON MULTI-PLATFORM GPS-ASSISTED TIME 
DIFFERENCE OF ARRIVAL SYSTEM 

A. PROJECT BASIS 

Given the problem of developing a portable, GPS assisted TDOA system with 
commercially available products, the project team at ARL:UT set out to define the 
problem and develop solutions based on which portions of the problem could be most 
influenced. The descriptions and performance reports below are based upon [15] and 
numerous conversations with team members. 

The TDOA problem can be reduced to three sub-problems: geometry, signal 
timing and measurement. All three are related to the rms error of the geolocation by 

Opasiuon « GDOP ■ & measurement (5-1 ) 

The Geometric Dilution of Precision (GDOP) is a scalar multiplicative factor 
which serves to magnify any timing and measurement error. It is directly related to the 
relative positions of the observers and the emitter. Poor GDOP(s), generally numbers 
greater than three in this application, result from geometries where the observers are lined 
up or clustered in one area with respect to the emitter. Conversely, good GDOP(s), 
numbers smaller than about three, result from situations where the observers are 
positioned such that their TDOA(s) intersect at nearly right angles. In three dimensions, 
the most favorable GDOP occurs by maximizing the volume of the polyhedron formed by 
pointing four vectors from the four observers to the emitter and enclosing the figure that 
results. 

The geometry of the situation is essentially out of the control of the TDOA system 
as the location of the emitter is not known a priori but would be needed in order to 
position the observers in an optimal configuration. Only when many more observers 
participate than the minimum needed for a determined solution (three are needed for 2D 
and four for 3D) can the geometry of the situation be better controlled by the TDOA 
system. In practice, this occurs rarely and thus, the best course of action considering (5-1) 
above is to minimize the timing and measurement error to the greatest extent possible. 
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GPS provides the capability to reduce both these error sources significantly. It was with 
this in mind that the ARL:UT project team, led by Mark Leach, pursued the design and 
implementation of a new GPS-assisited TDOA system. 

NSGSA in its tasking to ARL.UT added some additional requirements in the area 
of component sources, architectures and interfaces with existing equipment. The 
following sections outline the system in progressively more detail starting with a general 
summary and concept of operations, continuing into a hardware description and ending 
with a more detailed discussion of error sources and system performance in operational 
tests. 

B. PROJECT OVERVIEW 

r 

The system would receive only a frequency of interest in the HF, VHF or UHF 
line of sight communications frequency bands. It must utilize an open system 
architecture and maximize the use of commercial and government of the shelf hardware 
(COTS and GOTS respectively). It was permitted just a single voice grade channel for 
communication between observers and needed to interface with existing Navy hardware 
for that purpose (i.e. ANAVSC-3 transceiver). Finally, the system needed to be 
compatible with the Navy’s Unified Build (UB) environment and interact seamlessly with 
the Joint Maritime Command and Information System (JMCIS). 

The system in its current configuration consists of a network of observers 
(minimum of three for a geolocation in 2D and four for 3D) with one observer acting as 
the master node, the others as slave nodes. It employs the distributed processing 
technique depicted in Figure 5-1. Once given the frequency of interest by an existing 
signal acquisition system , the master node tunes its receiver to that frequency and 
notifies the slaves to do likewise. Each node samples the incoming signal and buffers the 
data in mass memory storage. The master logically determines the 400 ms portion of its 
samples that contains the most energy from the SOI and computes the Fast Fourier 
Transform (FFT) of this 400ms. The Master then broadcasts the FFT to the slaves with a 
GPS time stamp. 
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Figure 5-1 - Distributed processing, from [15] 

Once the slaves receive the FFT from the master they search their buffers for the 
400 ms of samples that correspond to the epoch of the time stamped FFT. They then take 
the FFT of their data and perform a CAF to extract the TDOA. In addition to this TDOA, 
the slaves send their GPS derived position, time and velocity (PVT) measurements 
corresponding to the time of intercept to the master. 

The master collects the TDOA(s) and PVT(s) and forwards the package to a 
standard Navy TAC3/4 workstation which calculates the geolocation. 

C. HARDWARE DESCRIPTION 

The system is comprised of two major sub-systems: the measurement subsystem 
and the geolocation subsystem. Each is designed in modular fashion accepting certain 
formatted inputs and producing certain formatted outputs. This modularity facilitates 
changes internal to either system including entire components or software modules as the 
process within each of these subsystems is inconsequential to the other module. Only the 
inputs and outputs are relevant. 

The measurement subsystem includes the hardware and software used to compute 
the TDOA(s). It is based on a 13 slot, C-size VXI chassis housing an HP Signal Analyzer 
package. This package includes a Motorola 68040 based VXI controller, a 10 MHz, 23- 
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bit (18 linear) analog to digital converter (A/D), a digital signal processing (DSP) card, 
200 MB of disk storage and software. The DSP module is a Motorola 96002 based 
floating point DSP chipset on a VXI card. Also on the chassis are a Ball-Efratom model 
FRK 10 MHz Rubidium (Rb) reference oscillator, a Precise Positioning Service (PPS) 
capable GOTS GPS receiver by Ashtech Inc. and E-systems and an ARL:UT developed 
clock synchronization module. The phase lock loop (PPL) module allows the Rb 
reference oscillator to be in synch with the GPS 1 pulse per second (pps) timing output. 
The bulk of this subsystem is shown in Figure 5-2. 




Figure 5-2 - Hardware configuration, from [15] 

The geolocation subsystem includes a TAC3/4 computer on which reside the 
algorithms both to compute the geolocation and communicate that fix to the JMCIS. In 
addition, it also holds the system’s capability for simulation and playback modes. The 
TDOA algorithms consist of the closed form solution presented in Chapter V and a 
Kalman filter method. The closed form or determined solution is the primary solution 
and lends itself quite naturally to situations where the emitter has short transmission 
durations which preclude the use of the Kalman filter method. In the event of long 
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duration transmission, the determined solution will feed the Kalman filter a first 
prediction and allow the Kalman filter to produce fixes as a function of time over the 
duration of the transmission. The Kalman filter, when given sufficient data and a starting 
point, has better error prediction capabilities and is thus the preferred method in the 
correct situations. 

The interface with JMCIS provides the system a data display capability. By 
updating the track database through a standard UB message, the system is able to both 
communicate with the Navy’s latest command and information system and avoid having 
to produce its own man-machine interface on this output level. 

Finally, this systemjs also capable of simulating scenarios in order to investigate 
the effects of variations on the many factors involved from the receipt of the signals to the 
processing of the data to the computation and display of the fix. In addition, a capability 
to playback sessions recorded and stored on tape is built in to the geolocation subsystem. 
This feature allows for post processing and analysis of operational tests enabling finer 
detailed analysis of many of the variables including error source investigation. 

Perhaps the one key point that stands above the other accomplishments of the 
system is the significant reduction of signal timing errors between the collection nodes. 
With the use of the PPS-capable GPS receivers, the 1 pps output from the GPS satellites 
in concert with differential GPS techniques enables the time offset between two GPS 
receivers, separated by up to hundreds of kilometers, to be on the order of 20 ns. In 
addition, the PPL circuitry allows each node’s Rb standard to be synchronized to within 5 
ns of this 1 pps output giving a total signal timing error of a mere 25 ns. In terms of light 
distance, that translates to 7.5 meters. Figure 5-3 shows the entire system’s functional 
relationships in summary. 
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Figure 5-3 - Tactical TDOA Functional Diagram, courtesy ARL:UT 



D. ERROR ANALYSIS 

Many sources of errors exist in the operational system. Another accomplishment 
of the project is the detailed identification and prediction of the magnitude of those error 
sources which have been validated by the operational tests outlined in the next section. In 
their investigations, the ARL:UT team was able to collect the sources into four major 
categories: fundamental limits, instrumentation limits, operational environment and 
network geometry. 

The lower bound on timing and measurement accuracy may be related to the 
variance of the TDOA estimator. This in turn is related to the incoming SNR, signal and 
noise bandwidth and integration time as first shown by Hertz and Azaria in [2] and later 
extended more specifically to this system by Rusu and Giulianelli in [16]. 

Instrumentation limits add significantly to the errors in signal arrival time 
determination. Specifically, though ideally each GPS receiver will output the GPS 1 pps 
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at the same time during each one second epoch, there will be some offset in reality. In 
addition, that 1 pps timing pulse is used to trigger the sampling of the incoming signal. 
Delays from the time the A/D converter receives the 1 pps and the onset of sampling for 
that clock cycle will vary from one A/D to another. This compounds the 1 pps offset 
error and translates directly to error in the cross-correlation computation of the TDOA(s). 
The cross-correlation methods of TDOA determination measure the time difference by 
determining a peak in the computed function. Five unknowns exists in this computation: 
signal frequency, signal bandwidth, noise bandwidth, SNR, and statistical properties of 
both signal and noise. The width of this peak and the presence of multiple peaks in 
indicative of the magnitude of the lower bound on the error in the measurement discussed 
above and is the direct result of a combination of these factors. 

The operational environment and network geometry constituted the bulk of the 
uncontrollable error contributions. While during testing, the environment is somewhat 
controlled, the operational system would certainly be subject to the conditions of the 
current tactical situation. This holds true as well for the geometries. 

E. PRELIMINARY TEST RESULTS 

Published test results include both static and dynamic emitters as targets. In 
November 1994, the geolocation of a static emitter in 2 dimensions was demonstrated in 
the Austin, Texas area using three observers with baselines of approximately 33 km. The 
target was a stationary NOAA FM transmitter at 162.4 MHz. The system at this point 
tested with a 2-D rms error of 124 m. Details can be seen in [15]. 

Dynamic tests were conducted in January 1995 using a moving emitter, two 
stationary observers and one moving observer. Again the tests were conducted in the 
Austin area and in just two dimensions only given just three observers. The target emitter 
was an aircraft at approximately 2000 feet, traveling at speeds of approximately 100 knots 
and transmitting with a 20 W UHF signal at 461.1 MHz. The moving observer was a van 
mounted node traveling at speeds of up to 35 mph in a 2 square mile area. During this 
test, the prototype system performed geolocations to within 277 m rms error using the 



49 



Kalman filter solution. This accuracy was achieved despite various geometries that 
ranged from outstanding GDOP(s) of almost 1 to very poor GDOP(s) of 10. 

While these tests failed to achieve the goal of 100 m rms error, instrumentation 
changes promised to increase the accuracy without significant changes to the system. 

Two additional tests, each more extensive and demanding than the previous, have been 
conducted with the system. While no performance results have been formally published, 
shipboard tests off the coast of San Diego California have proven the system capable of 
120 m rms error against representative VHF/UHF dynamic emitters with excellent 
geometry and ~400 m rms error with poorer geometries [17]. Further changes in 
receivers promise to remove biases in the measurements created by the current radios and 
improve performance proportionally. 

In addition, the use of a cyclostationary TDOA determination algorithm should 
provide some improvement in the performance of the system low SNR environments. 

The initial step for the incorporation of cyclostationary signal processing into the 
ARL:UT system is to code the algorithms in MATLAB® and test them with ARL:UT test 
data. The coding of the algorithms in MATLAB® appears in Chapter VI. 
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VI. ALGORITHMS 



A. BACKGROUND 

Any implementation of a TDOA processor utilizing SPECCOA must provide both 
an SCF as well as a CSCF at a cycle frequency oto, characteristic of the SOL In [19], 
Loomis and Berstein develop a functional TDOA processor using SPECCOA that 
provides for the efficient computation of both spectral correlation functions as well as the 
selection (Xo- This system represents the basis for the implementation of SPECCOA and 
the geolocation algorithm in MATLAB® for this project and appears in Figure 6-1. 




Figure 6-1 - TDOA geolocation processor functional diagram, adapted from [19] 

The four blocks of the TDOA processor constitute progressively lesser degrees of 
computational complexity. The most complex of the operations is the calculation of the 
SCF(s) in the first operation of the processor. The SCF(s) of the master and the four 
slaves are necessary for the selection of the appropriate cycle frequency in the next 
operation of the processor, the characterization. The SSCA algorithm computes the four 
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SCF(s) and uses a plotting routine to display them. In the next operation, the 
characterization is performed manually using a plot of the maximum values of the SCF(s) 
for each value of cycle frequency. This plot facilitates the choice of the optimum cycle 
frequency for use in SPECCOA. The third operation involves the computation of the 
SCF and CSCF at the cycle frequency selected during characterization. Computed with a 
simple frequency smoothing algorithm, the SCF and CSCF at the cycle frequency of 
interest feed the final TDOA processor operation, SPECCOA. SPECCOA in turn 
determines the TDOA by plotting the solution to the ad hoc MLS optimization found in 
(3-32). The maximum value of this function represents the TDOA. 

The SCF(s) in the first operation as well as the characterization must be computed 
for all four observers for each set of data. Cycle frequencies for the same emitter vary 
with observer due to instrumentation mismatches at each location. Details of cycle 
frequency selection in these situations appear in Chapters VII and VIII. The cycle 
frequency specific SCF and CSCF are computed just three times because as noted in 
Chapter V, the only TDOA(s) used in practice are those between the master node and the 
three slaves. No cross slave TDOA(s) are computed or used in the prototype system and 
thus they are not computed here. Finally, three runs of SPECCOA are required to feed 
the geolocation algorithm. 

B. STRIP SPECTRAL CORRELATION ANALYZER IN MATLAB® 

The SCF(s) computed in the first operation of the TDOA processor can place a 
significant burden on the processor if not optimized for computational efficiency. In 
order to achieve the efficiency necessary for use in a tactical environment, a highly 
expeditious method for computing and plotting the SCF is used. Directly from its 
definition, (3-17), the SCF can be computed using a traditional time-smoothing approach 
as depicted in Figure 6-2. While this approach provides the most accurate estimate of the 
SCF, it proves prohibitively time consuming and inappropriate for most applications [20]. 
The SSCA algorithm, a modification of the time-smoothing method, provides the best 
combination of accurate estimation of cycle frequencies of interest and computational 
efficiency [21]. By eliminating one of the band-pass filters in Figure 6-2 and using a 
Fourier transformer at the input and output, SSCA computes the SCF along 
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Figure 6-2 - Time-smoothing spectral correlation analyzer [20] 

diagonal strips of the bi-frequency plane as depicted in Figure 6-3. While the SSCA does 
manage to compute the SCF in the most efficient manner developed thus far, it has some 
compromise. The output signal to noise ratio suffers slightly as a result of the efficiency 
gained by the elimination of the filter and the use of the Fourier transform. However, the 
computational savings gained far outweigh the minor degradation in signal power to 
noise power at the output of the analyzer. 

Figure 6-4 illustrates the SSCA architecture. The input filter consists of a sliding, 
Hamming-windowed coarse FFT of length N ' . Each of these N' length segments is 
essentially the output of a band-pass filter with a bandwidth of 1/ N ' . This band-pass 
output is then downconverted to baseband by a complex exponential multiplication. 
These products are termed the complex demodulates of the input signal and appear in 
greater 
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Figure 6-3 - SSCA computation in the bi-frequency plane, N'= 8 




Figure 6-4 - SSCA architecture, from [20] 
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detail in Figure 6-5. As can be seen, the input sequence is essentially decimated by the 
process of windowing and sliding at the input filter. Loomis and Brown found that the 
most appropriate value for the decimation factor or subsampling parameter L was N' / 4 
[ 20 ], 



X(n) (length N) 




O 



o 



N' 



points each (Hamming windowed and FFT’d) 



o 



o 



o 

origin slides by L each time 



N/L 



Complex Demodulates 

- products of above segments 
and complex exponentials 



N 




Figure 6-5 - Detail of the computation of complex demodulates 

Because these complex demodulates are multiplied by the original sequence 
before the final FFT, they must be interpolated to the original sampling rate. In [20] and 
[21] Loomis and Brown use a simple hold operation. Linear interpolation induces less 
high frequency error and adds little in computational complexity. It is used here in place 
of the first order sample and hold interpolation. Finally, the interpolated complex 
demodulates and the original sequence are complex multiplied and the product Fourier 
transformed by a full length, unwindowed FFT. The relationship of the output to the bi- 
frequency plane appears in Figure 6-6. Each row of output represents one of 
the N' diagonal strips in the bi-frequency plane as depicted on the previous page in Figure 
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6-3. Each column represents the values of the SCF along the diagonal lines. The 
resolution along lines of /is coarse at 1 / N' as a byproduct of the computational savings. 
However, the resolution in a, at 1/N, is much higher which accommodates the 
characterization nicely. 




Figure 6-6 - Block form of SSCA output 

Given that the majority of operations in the SSCA can be implemented as vector 
or matrix operations, its implementation in MATLAB® is fairly simple. The psuedocode 
for the algorithm appears in Figure 6-7. The MATLAB® code for the computation of the 
SCA over the entire bi-frequency plane, called ssca, appears in Appendix A. Of note, 
several MATLAB® functions forced the use of additional statements that might not 
otherwise be found given the psuedocode. Th ejft function in MATLAB® computes both 
positive and negative frequencies but does not place the origin in the center of the output 
vector. Instead, the output of fft takes advantage of the cyclic nature of the DFT. Its 
output begins at the origin of the frequency, moves through f s /2 and then mirrors the 
negative frequencies with the positive continuing at -f s /2 and ending with zero again. 
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This forces any application needing a double sided spectrum to rearrange the vector so 
that the negative frequencies start the output rather than end it. MATLAB® has an 
fft shift function that performs such an operation and appears in the code in several places 
for that reason. The output of ssca is a matrix with N' rows corresponding to 
the N' diagonal strips of the computation and columns corresponding to the SCF along 
the respective strips. Referring again to Figures 6-3 and 6-5, these strips are not 
computed along lines of / and a and must therefore be transposed along those lines before, 
plotting. 



/* Compute Complex Demodulates of x */ 

Do p:= 0 to P-1 

Compute x T (pL,f k ) = FFT[a(r)x(pL+r)] 

Do k := —N'/2 to AC/2-1 
Compute X T (pL, f k ) = xT(pL, f k )exp{-j27tpLk/N' } 

end 

end 

/* Interpolate X T (pL, f k ) */ 

Compute X T (pL, f k ) => X T (n,f k ) 

/* Compute and Smooth Product Waveforms 
Dok := -N'/2to N'/ 2-1 
Compute S£ r (n,f k ) = y* (n)X T (n,f k ) 

Compute S“; (n, /, )* = FFT{g(n)S^ {n,f„ )} 

end 



Figure 6-7 - Psuedocode for SSCA computation, from [21] 

The program plotscf performs that function transposing the matrix output of the 
ssca program and plotting both the 3-D SCF and the 2-D magnitude of SCF versus a for 
the characterization operation. The transformation involves applying the relationship 
between the diagonal strip subscript k, the frequency/ and the cyclic frequency a. Those 
relationships from [22] are 
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( 6 - 1 ) 



fk= k 




2 2 



a = 2f t -2f 



( 6 - 2 ) 



Given these relations and the MATLAB® s«//function, the SCF plots easily 
along lines off and a. The 2-D plot for additional characterization aid follows directly 
from the SCF plot. 



C. SPECTRAL COHERENCE ALIGNMENT IN MATLAB® 



SPECCOA represents yet another algorithm highly suited for MATLAB®. The 
majority of functions associated with computing (3-32) are vector operations which 
simplifies the MATLAB code considerably. Figure 6-8 is a basic block diagram of 
SPECCOA. 




Figure 6-8 - SPECCOA block diagram 
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Once the characterization has determined the cycle frequency of interest, oto, a 
frequency smoothing algorithm estimates the SCF and CSCF for oto. The inverse Fourier 
transform of the complex product is then multiplied by the kernel e ,Ka * x . The maximum 
value of the real portion of this product represents the TDOA of the SOI between the 
observers. 

SPECCOA is used three times to determine the TDOA between the master node 
and the three slave nodes. These TDOA(s) along with observer positions are inputs to 
Dr. Rusu’s geolocation algorithm the output of which is the location of the emitter. 

D. THE TDOA CLOSED FORM SOLUTION IN MATLAB® 

Dr. Rusu’s closed form solution lends itself readily to MATLAB implementation. 
The inputs include four observer position vectors in WGS84 coordinates and three 
TDOA(s), D 12 , D 13 and D 14 . All the intermediate variables are either vectors of length 3 
or matrices no larger than 3x3. The final solution involves finding the roots of a 
quadratic that represent the time of emission of the SOI and substituting that time into a 
range equation to find the WGS84 coordinates of the transmitter. 

Of particular note to the MATLAB® code that appears in Appendix A is the use 
of the backslash function, \, in lieu of the matrix inverse in the algorithm. This function 
was specifically developed for use in solving equations of the form 



Ag = q (6-3) 

for g where A is a matrix and g and q vectors. Normally, in matrix algebra the solution 
has the form 



g = A-'q 



(6-4) 



if A has an inverse. However in MATLAB®, it is more accurate to use the matrix 
division function or backslash. The accuracy of the backslash in this scenario approaches 
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machine precision while the use of the matrix inverse can reduce that accuracy several 
orders of magnitude [23]. 

One problem persists with the algorithm but is a function of the closed form 
solution rather than its implementation. The two roots both represent viable solutions to 
the problem. Information outside the algorithm must be used to eliminate the ambiguity. 
Multiple solutions from moving observers would more often than not cause one solution 
to converge while the other diverged. No matter how its resolved outside the algorithm, 
the ambiguity always poses a problem insurmountable by the mathematics. 

The algorithms need to be tested on controlled data prior to evaluation with 
ARL:UT test data. Several signals were generated using Simulink®. These signals serve 
to ensure the proper functioning of all the code written for this work. Following the 
generated data, ARL:UT data provides a vehicle for comparing the current CAF based 
ARL:UT method with SPECCOA. The test plan is detailed in Chapter VII. 
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VII. TEST PLAN 



A. GENERATED SIGNALS 

Three signals are used in initially testing the TDOA processor. Generated with 
Simulink®, all three are BPSK signals and vary from noiseless to very poor SNR. Each 
signal was sampled at 100 kHz. They had carrier frequencies of 0.25 f s and data rates of 
0.05 f s . All were modulated with a random binary signal generator. They are used in 
various combinations to verify the functions of the three major programs, ssca, plotscf 
and speccoa. 

The first set of tests involves using a noiseless BPSK signal to first verify ssca 
and plotscf. Given that these programs function properly, speccoa is tested finally. 
Figure 7-1 shows the simulation used to generate noiseless BPSK. A total of 1 second of 
signal is actually collected and stored as testbpsk.mat. 









Random Binary 
Wave 




BPSK Modulation 




testbpsk.mat 













Figure 7- 1 - Generation of noiseless BPSK as testbpsk.mat 
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This signal approaches the ideal for the processor for several reasons. First, it is 
completely binary from its inception. No analog to digital conversion is necessary and 
thus introduces no error in the process. Second, it is completely noiseless. Finally, it is 
modulated by a random binary wave which introduces no biases in the modulated 
waveform. Given this, its SCF should approach that given by theory taking into account 
the time smoothing effects of the SSCA algorithm. Figures 7-2 and 7-3 are the time 
domain and frequency domain plots respectively for testbpsk.mat. 




time (s) 

Figure 7-2 - Time domain plot of testbpsk.mat 

The ideal SCF of a noiseless BPSK signal can be seen in Chapter III, Figure 3-4. 
The data feature of testbpsk.mat should appear at 0.05 / s while the twice carrier feature 
in conjunction with the data feature appear at 0.45/ s , 0.5/ s and 0.55 / s . Figure 7-4 
displays the results of ssca and plotscf for N = 4096 point sample. Only the positive 
cycle frequency portion of the bi-frequency plane is plotted. The cycle frequency axis 
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draws from upper right to lower left and is plotted in terms of nf s /N, n as axis label, 1 IN as 
cycle frequency resolution. The spectral frequency follows from upper left to lower right 
labeled in the same manner. They are unlabeled in the figure to eliminate the additional 
clutter. As noted in Chapter VI, the SSCA algorithm trades resolution in the frequency 
domain for computational efficiency without losing the four key features. 
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Figure 7-3 - Frequency domain plot for testbpsk.mat 

Finally for testbpsk.mat, speccoa is run with observer 1 collecting testbpsk.mat 
samples 1 to 8192 and observer 2 collecting testbpsk.mat samples 71 to 8263 for a 70 
sample delay. The results of this test appear in Figure 7-5. Because the signals at the two 
observers are completely correlated, this scenario represents an ideal situation for the 
algorithm. Its complete success is expected for this case but serves as an initial check of 
speccoa. The remaining two signals introduce both noise and poorly correlated signals to 
the scenario. 
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Figure 7-4 - SCF (top) and Mag(SCF) vs alpha (bottom) for testbpsk.mat 
computed and plotted by ssca and plotscf 
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Two very noisy BPSK signals constitute the remaining generated signal testing for 
the system. They were created by passing the noiseless BPSK signal through separate 
AWGN channels with Eb/No of 0.88 or -0.534 dB. This corresponds to probability of bit 
error of 10 2 . Figure 7-6 shows the generation of nlbpsk.mat and n2bpsk.mat. Because 
the AWGN channel introduces significant noise the BPSK signals, nlbpsk.mat and 
n2bpsk.mat are highly corrupted. Figures 7-7 and 7-8 show the time domain and 
frequency domain plots respectively for the two signals. 
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Figure 7-5 - Output of speccoa routine for testbpsk.mat with 70 sample delay 

Naturally, the level of noise in these BPSK signals will be evident in their SCF(s). 
Figures 7-9 and 7-10 display the output of ssca and plotscf. Note that the reduction in 
SNR compared to testbpsk has added noise to the plots but has not masked the important 
features. While in theory as presented in Chapter III, the AWGN should have no features 
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other than those evident at a =0, it does add noise to the higher values of a in the SCF 
plots in Figure 7-9. It does so because the AWGN channel used in the simulation is not 
ideal AWGN and thus is not truly random and Gaussian. Thus, some statistical 
impurities in the channel as well as the computational trade-offs in the estimation 
procedure of the SSCA algorithm lead to contributions above a = 0. 




Figure 7-6 - Generation of corrupted BPSK as nlbpsk.mat and n2bpsk.mat 

Finally, speccoa uses nlbpsk.mat for observer 1 and n2bpsk.mat delayed by 70 
samples for observer 2 to verify that it indeed works in highly corrupt environments given 
the proper signal features as inputs. The results appear in Figure 7-11. The successful 
conclusion of testing with generated signals leads next to testing of the ARL:UT data 
from a November 1995 test in the Austin, Texas area. 
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Figure 7-7 - Time domain plots for nlbpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-8 - Frequency domain of nlbpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-9 - SCF(s) for nlbpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-10 - Mag(SCF) vs. alpha for nlbpsk.mat (top) and n2bpsk.mat (bottom) 
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Figure 7-1 1 - Output of speccoa routine for nlbpsk.mat and n2bpsk.mat 

B. ARLrUT COLLECTED DATA 

The portion of the November 1995 Austin tests that are used in this work were 
collected by four stationary observers in the Austin area. The emitter was a push-to-talk 
narrow-band FM (NBFM), 1 W transmitter located near the center of the four observers. 
The GDOP of the scenario, given the positions of observation and the stationary target, 
was excellent at less then two. 

During the test, each observer records -800 ms of data every 5 seconds. This data 
was stored to magnetic tape for post test analysis and further off-line testing. It is used to 
compare both TDOA determination and geolocation performance of the ARL:UT 
algorithms and those developed here. In order to best understand the output of ssca, 
plotscf and speccoa for the ARL:UT data, it is necessary to follow the signals’ processing 
from antenna to storage media. The processing at each observer is identical with respect 
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to archiving the data for later use and thus a single explanation serves for all four 
observers in this case. 
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Figure 7-12 - Beginning of ARL:UT signal processing 

Following reception, the 10 kHz bandwidth signal is heterodyned to an 
intermediate frequency (IF) of 21.4 MHz. At the IF, the 10 Msamples/second A/D 
converter aliases the signal at +1.4 MHz and -1.4 MHz with additional spectra added at 
intervals of the sampling rate. The spectrum at 1.4 MHz is downconverted to baseband 
as the product of an appropriate complex exponential multiplication. Its mirror at -1.5 
MHz is shifted to -2.8 MHz in the process. This series of steps is summarized in Figure 
7-12 with a spectrum relatively representative of that of the test data. 

At this point, the data is decimated in time by a factor or 1024. This step again aliases the 
signal at intervals of the new sampling rate,/ s = 9765.625 Hz. The result is a set of four 
spectra in the interval -/ s /2 to + / s /2. Figure 7-13 shows the results of this process in the 
time domain and frequency domain. Note the relationship between the highest negative 
frequency spectral line at approximately -440 Hz and the highest positive frequency 
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spectral line at approximately 4440 is exactly /s/2 as they are related, aliased versions of 
the spectrum. The lowest negative frequency at approximately -4460 Hz and lowest 
positive frequency at approximately 420 Hz share that same relationship for the same 
reason. The frequencies are approximated above as they change slightly throughout the 
test and vary by observer as well. 

The corresponding SCF has features along the cyclic frequency axis that 
correspond to the various combinations of separation between the four spectral lines in 
the PSD of the signal. The greatest contribution to the SCF occurs at a separation of/ s /2 
given the two constant relationships between the spectral lines sited above. At a =/ s 12, 
all four lines contribute to the SCF thus the a =/ s /2 feature contains the most energy. 

This is best illustrated the SCF and the magnitude vs. cycle frequency plots in Figure 7- 
14. Note the greatest peak in the SCF occurs at/ s /2. 

Finally, using this greatest feature, speccoa finds the TDOA between observers. 
While other features may also lead to reliable answers, the f s H feature is a constant 
throughout the testing. All observers have this feature and it is always the strongest 
feature of the SCF. An example of the output from speccoa appears in Figure 7-15. It 
depicts the TDOA between the master node and slave #1. Specific results from the testing 
including some additional steps taken for more resolution in the TDOA computations 
appear in Chapter IX. 
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Figure 7-13 - Time (top) and frequency (bottom) plots of master signal, time 409400 
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Figure 7-14 - SCF (top) and Mag(SCF) vs. alpha for master signal, time 409400 
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Figure 7-15 - Output of speccoa for master and slave #3, time 409510 
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VIII. GEOLOCATION SOFTWARE WORKBENCH BLOCKS 



A. SIMULINK® BLOCK CONSTRUCTION 

In addition to the testing of the combination of SPECCOA and Rusu’s closed 
form TDOA geolocation solution, the algorithms serve as the first blocks readied for 
construction for a developing geolocation software workbench. The details of their 
construction are outlined below in theory. 

Blocks that may be used in simulations run in the Simulink® environment 
conform to a general, highly flexible standard creates by The Mathworks, Inc. 
programmers. They may be implemented in one of three ways: graphically, through ni- 
fties or through another language such as ANSI C which produces a mex-file. No matter 
the method, the result is an S-function that operates according to the rules and protocols 
of the Simulink® environment. [24] 

Graphical creation of an S-function involves connecting the existing blocks that 
perform the needed function, collecting them as a single entity and masking that 
collection as the new block. This procedure forces Simulink® to create an S-function 
that describes the new system. This is by far the easiest method of creation for a speccoa 
block and a closed form TDOA geolocation block. The other two methods, m-files and 
mex-files, are appropriate in some instances such as certain state-space systems or other 
applications where the system can easily be written as a set of differential equations. 
Though more complex to create, S-functions created as m-files and mex-files simulate 
considerably faster. However, with the addition of the Simulink Accelerator®, 
graphically created S-functions enjoy some improvement in execution time as well. [24] 

B. SPECCOA BLOCK 

All the functions needed to use the speccoa code written for the testing described 
in both Chapters VII and IX, exist in Simulink® or one of the toolboxes. The Signal 
Processing Toolbox ® and DSP Blockset® contain the necessary buffers and functions 
needed by the speccoa code. A user defined Matlab® function block serves as the 
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vehicle for the speccoa code to be used in the Simulink® environment. Connecting the 
buffers with the user defined function block constitutes the majority of the speccoa block 
In general, a geolocation workbench block of speccoa requires two scalar 
sequences, the time samples of the two signals and one scalar parameter, the cycle 
frequency of interest as inputs. Time domain samples feed a buffer of length defined by 
the user in order to achieve the desired resolution. These time domain samples come to 
the speccoa block from an appropriately designed data conversion block. Regardless of 
its input, the conversion block must provide the speccoa block with single real or 
complex values for each simulation time period. In addition, the characteristic cycle 
frequency enters the block following the operation of a characterizer (defined elsewhere) 
to complete the inputs. The output of the speccoa block is a 2 x 1 vector, a time-stamp 
and a TDOA. Figure 8- 1 shows a schematic of the speccoa block. Its role in a possible 
simulation scenario appears in Figure 8-3. 




Figure 8-1 - Schematic of speccoa geolocation workbench block 
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c. 



GEOLOCATION BLOCK 



The closed form TDOA geolocation block is the more complex of the two 
algorithm blocks. Its inputs include four 3 x 1 observer position vectors with x, y and z 
parameters for each observer. In addition, it must have the output of three speccoa 
blocks. The geolocation block’s output is a 4 x 1 vector that includes the target’s x, y and 
z coordinates and a time corresponding to that geolocation. Figure 8-2 illustrates a basic 
diagram of the block. Its role in a possible simulation scenario appears in Figure 8-3 as 
well. 




Figure 8-2 - Closed form TDOA geolocation solution block schematics 



1). SIMULATION SCENARIO 

One possible scenario that might be simulated using the two blocks described 
above is the November 1995 Austin test of the ARL:UT system. Using the signal data 
files and position files for each of the four observers, both the speccoa block and the 
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closed TDOA form geolocation block are testable. The addition of a characterizer is 
necessary to feed the speccoa block. This function could be performed by a user defined 
block that calculates the maximum value of the SCF for each value of cycle frequency as 
plotscf plotted (see Chapter VII). The maximum value of this function could be the 
characteristic cycle frequency used in the speccoa block. In addition, the geolocation 
block could feed a plotting routine to display each geolocation. A diagram of this 
implementation appears in Figure 8-3. 




Figure 8-3 - Simulation scenario 

Clearly the Simulink® environment is ideal for testing the many aspects of 
geolocation from data processing to parameter computation to geolocation estimation. 
The user defined MATLAB® function block provides for easy introduction of m-file 
routines into the S-function domain. It allows the interchange of various blocks in any 
scenario. This fact alone lends great flexibility and power to a geolocation software 
workbench centered around this environment. Other blocks such as pre-filters for the 
data or logic to evaluate the quality of the TDOA(s) or geolocation can also be 
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constructed. In addition, data from many different sources is handled simply by 
modifying the data conversion block to accept a different form of input. Its output 
remains the same. Using the Simulink® environment eliminates the need to define 
standards for module construction and provides an established simulation engine with 
many options. It is surely the most expeditious and intelligible implementation of a 
geolocation simulation system. 
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IX. RESULTS 



A. BACKGROUND 

Two comparisons are made between the cyclostationary TDOA processor 
measurements and those of the ARL:UT system. Both the TDOA(s) and the resulting 
geolocations serve as benchmarks of performance. In comparing the two systems’ 
results, theoretical TDOA(s) calculated from the target’s position and the positions of the 
four observers, using a value of 2.998 x 10 s for c, the speed of propagation, are the 
baseline for comparison of the TDOA measurements. The target’s actual position 
throughout the test acts as the baseline for comparison of the geolocations. 

The TDOA(s) measured by both systems are biased by the radios used in the test. 
This constant bias was discovered in post test analysis and removed from the ARL:UT 
published results. However, the biases remain intact here to simplify the calculations 
and comparison. The geolocations suffer from the same biases and are also left 
uncorrected for the comparison. 

The comparison consists of 25 time periods of data from the November 1 995 
Austin test. Because all four observers recorded data during these times, geolocations in 
three dimensions are possible. The process of arriving at three TDOA(s) and a 
geolocation involves plotting the SCF(s) for the four observers, choosing a cycle 
frequency, feeding that data to speccoa and finally running the closed form solution for 
the geolocation. The plots for the four SCF(s), the four magnitude of SCF vs cycle 
frequency and the three speccoa outputs appear in Appendix B for each of the 25 time 
periods. As noted in Chapter VII, the feature of the SCF at fjl provides speccoa the 
needed cycle frequency in every time instance. A summary of the results appears below. 

B. TDOA COMPARISON 

The comparison of TDOA(s) between the master node and slaves 1 , 2 and 3 
appear below in Figures 9-1 through 9-6 below. Figures 9-1, 9-2 and 9-3 merely plot the 
TDOA(s) in the order in which they were determined. 
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Figure 9-1 - TDOA(s) for master node and slave #1, ARL:UT (top), speccoa (bottom) 
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Figure 9-2 - TDOA(s) for master node and slave #2, ARL:UT (top), speccoa (bottom) 
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Figure 9-3 - TDOA(s) for master node and slave #3 ARL:UT (top), speccoa (bottom) 
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Figures 9-4, 9-5 and 9-6 show the deviation in seconds, from the calculated, 
theoretical TDOA, for each of the 25 time periods along the horizontal axis. The top 
plots in each figure again show results from the ARL:UT test; the lower plots show 
results of the cyclostationary TDOA processor. 

C. GEOLOCATION COMPARISON 

While the algorithm in both instances is identical, a comparison of the results 
serves to further illustrate the performance of the two different TDOA determination 
algorithms. The method of comparison is merely the geolocation performance. 
Unfortunately, even given the wide variety of TDOA(s) determined by the 
cyclostationary TDOA processor, the version of the closed form TDOA geolocation 
solution coded in MATLAB® produced nearly the same result for all 25 time instances. 
While this consistency is commendable, the various TDOA(s) should have produced 
fixes separated by several kilometers. This lack of definite variation in the geolocations 
reflects poorly on the quality of the algorithm’s outputs. 

The actual location of the transmitter is given by the coordinates 
x, = -741941.474 m, y, = -5462120.413 m and z, = 3198242.718 m. The MATLAB® 
solution produces xm = -740591 m, yM = -5443692 m and zm = 3187478 m. The only 
change to these rounded integer values over the course of the 25 time periods occurs in 
the decimal places. Clearly this is a very poor fix in addition to the lack of variation. 

The ARL:UT system successfully found the coordinates to within 100 m on all three axis 
for all 25 time periods. More often than not, the errors in the ARL:UT calculated values 
were on the order of 10 to 20 meters. 

Because the MATLAB® algorithm did not perform properly, a valuable 
comparison is impossible. The TDOA(s) discussed above are the only viable measure of 
comparison available. 



87 



deviation (seconds) deviation (seconds) 



xIO 



-7 




10 15 

times 



8 
6 
4 
2 
0 
-2 
-4 

0 5 10 15 20 25 

times 

Figure 9-4 - Deviation from theoretical for TDOA between master node and slave #1, 
ARL:UT (top), speccoa (bottom) 
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Figure 9-5 - Deviation from theoretical for TDOA between master node and slave #2, 
ARL:UT (top), speccoa (bottom) 
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Figure 9-6 - Deviation from theoretical for TDOA between master node and slave #3, 
ARL:UT (top), speccoa (bottom) 
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D. COMMENTS 



Clearly the cyclostationary TDOA processor does not perform as predicted by 
theory and as verified by simulation. Due to the decimation of the original time domain 
signals, several additional processing steps are necessary in order to achieve the needed 
resolution in the TDOA computation. These steps include resampling the signal at a 
much higher rate in order to measure a TDOA on the order of 10' 6 . This step causes 
problems for the speccoa code. 

The average values for the two TDOA determination methods appear in Tables 
9-1 to 9-3 below. Also in the table are the values for the average deviation of the 
TDOA(s) from the theoretical values. 

The code for the geolocation algorithm obviously has flaws as well. Its inability 
to distinguish between the widely varying sets of TDOA(s) with widely varying 
geolocations indicates a fundamental programming flaw. 



Algorithm 


TDOAmi 


Theoretical TDOAmi 


Deviation M i 


Cyclostationary TDOA 


-1.308e-5 


1.434e-6 


1 . 1 65e-5 


ARL:UT System 


-1.805e-6 


1.434e-6 


3.66e-7 



Table 9-1 - TDOAmi summary 



Algorithm 


TDOAm 2 


Theoretical TDOA M 2 


Deviation M 2 


Cyclostationary TDOA 


8.44e-7 


3.828e-6 


2.984e-6 


ARL:UT System 


2.697e-6 


3.828e-6 


1.13 le-6 



Table 9-2 - TDOAm 2 summary 
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Algorithm 


TDOAm3 


Theoretical TDOA M 3 


Deviation M3 


Cyclostationary TDOA 


-5.869e-6 


-3.124e-6 


2.745e-6 


ARL:UT System 


-3.47 le-6 


3.124e-6 


3 Alt-1 



Table 9-3 - TD0A M 3 summary 
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X. CONCLUSIONS 



A. OBSERVATIONS 

The cyclostationary TDOA processor in combination with the closed form 
geolocation solution in MATLAB® does not outperform the CAF based system. 

Because the results deviate significantly from the theory, several additional tests are 
necessary to isolate the cause of the anomaly. First, speccoa must run with fully 
decimated ARL:UT data. Because the decimation exceeds the resolution needed to 
determine the TDOA(s), properly the result should be zero in every instance. It is not. 
The algorithms returns random, obviously unreliable answers for the real world test data. 

The next step is the verification of zero delay with generated data. The 
generated BPSK signals with zero delay produce zero as a result with speccoa. Thus, 
speccoa indeed functions properly even at delays below the resolution of the sampling 
rate. However, it does not function properly with the ARL:UT data. 

Finally, since the SPECCOA algorithm reduces to a GCC at a = 0 in the absence 
of Doppler shift, speccoa can be tested as a GCC to reproduce the ARL:UT results. This 
test produces random, unreliable results yet again with ARL:UT data. It does, however, 
produce the correct results with the generated BPSK signals once again. 

To summarize, plotscf and the cycle frequencies it produces are verified by the 
correct plots of the generated signals and the correct determination of the dominant 
cyclic features of the BPSK signals even in the presence of significant noise. Thus, the 
cycle frequencies used in speccoa for both the generated signals and the ARL:UT data 
are deemed reliable. Speccoa on the other hand, produces perfect results for the 
generated signals at all significant cycle frequencies and at a = 0. However, it does not 
produce consistent results for the dominant cyclic feature in the ARL:UT data; nor does 
it produce results at other cycle frequencies that are consistent with the /s/2 feature’s 
results. 

These additional details lead to one most likely conclusion. The ARL:UT data 
does not reach the algorithms in the correct form from storage to processing. Most 
probably, the format is incorrect. In its reading into the MATLAB® environment, the 
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data is essentially corrupted producing such seemingly contrary results. The generated 
signals are generated in the MATLAB® environment, producing no such format 
problems. The ARL:UT data is originally stored in UNIX format, byte reversed for PC 
usage and then loaded into the MATLAB® environment for processing. It is this last 
step that is believed to cause the anomalies in the results. 

As a positive outcome, the testing did provide an idea for the computational 
complexity that a cyclostationary TDOA system would require. The process from 
loading the signals to finding a geolocation required approximately 2 minutes per 
geolocation on a Pentium Pro™ based machine with 64 MB EDO RAM. This is not 
tactically feasible. While coding the system in C would increase the speed significantly, 
the computations of the system would remain the limiting factor in its implementation. 

The cyclostationary TDO£ processor holds promise for near future 
implementation. The problems encountered are tractable. All that remains is the 
computational power to ran the system in a tactical environment. 

B. AREAS OF FURTHER STUDY 

Determining the proper method with which to load the ARL:UT data into the 
MATLAB® environment remains the highest priority for solving the major problems of 
the current algorithms. If this is the cause of the significant inconsistency encountered, 
its solution will surely improve the results to accuracies predicted in theory. 

The polynomial fit used by the resample function in MATLAB® needs to be 
examined closely to determine the extent of phase error it introduces to the resampled 
signal. This would allow for the development of a more accurate method of 
interpolating back to a sampling rate that provides sufficient resolution for the TDOA 
measurement. 

The closed form solution error must be found and corrected. It is likely that the 
error involves the reintroduction of c = 2.998 x 10 8 once the contributions of the 
TDOA(s) are mapped to the three coordinates. In his solution, Rusu sets c = 1 for 
simplicity. While this works for the derivation, it most likely must be given its proper 
value before the final solution can be reached. 
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Once the two algorithms were corrected, using real world signals from a more 
hostile signal environment for comparison would highlight the benefits of 
cyclostationary processing best. 

Finally, further development of Geolocation Software Workbench blocks would 
accelerate the introduction of a simplified testing platform for all aspects of geolocation. 
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APPENDIX A. MATLAB CODE FOR ALGORITHMS 



The programs that follow constitute the major sections of code used in this work 
They all contain headers and commenting for explanation purposes. They are described 
briefly in Table A-l below. 



Program Title 


Purpose 


loaddata.m 


Loads headers and signal data, converts 4 bytes real, 4 
bytes imaginary to complex vector notation. 


ssca.m 


Strip-spectral correlation analyzer. Calculates the SCF 
of the input sequence along diagonal lines according to 
the SSCA algorithm. 


plotscf.m 


Transforms the diagonal lines of ssca.m into lines along 
f and a. Plots both the SCF and the magnitude of SCF 
vs cycle frequency. 


speccoa.m 


Resamples time domain data to user specified value (p). 
Calculates the TDOA of two input sequences in terms 
of samples using SPECCOA algorithm. 


rusugeo.m 


Takes positions of observers and three TDOA(s) to 
produce geolocation in WGS84 ECEF coordinates. 
Finds both solutions. User must eliminate one. 




Table A-l - Summary of code 
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<^****************************************:f::************************* 

% 

% Name: David A Streight 

% 

% Naval Postgraduate School - Monterey, California 

% 

% loaddata.m 

% 

% Operating Environment 
% Windows 95 

% 

% Description 

% Loads 32-bit float ARL:UT data into four header and data vectors. 

% Master site is obs 1 , other sites are 2, 3 and 4. 

% 

% Date of last revision 
% 20 December 1 996 

% 

% Inputs 

% Four ARL:UT stored files 

% 

% Outputs 
% None 

% 

clear all; 

% open four observer files read only 
obsl = fopen('F:\data\nov_ut~2\mcc\409400.rev','r'); 
obs2 = fopen('F:\data\nov_ut~2\berg\409400.rev’,'r'); 
obs3 = fopen('F:\data\nov_ut~2\beecaves\409400.rev','r'); 
obs4 = fopen('F:\data\nov_ut~2\expo\409400.rev','r'); 

% load 3 header variables into header vectors 
[headerl] = fread(obsl,3,'float'); 

[header2] = fread(obs2,3,'float'); 

[header3] = fread(obs3,3,'float'); 

[header4] = fread(obs4,3,’float'); 

% load data into 4 data matrices col#l - real col#2 - imaginary 
[xl] = fread(obs 1 ,[8 192,2],'float'); 

[x2] = fread(obs2,[8192,2],'float’); 

[x3] = fread(obs3,[8192,2],'float'); 

[x4] = fread(obs4,[8 1 92,2], 'float'); 
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% convert data matrices to complex vectors 
si =xl(:,l)+j*xl(:,2); 
s2 = x2(:,l)+j*x2(:,2); 
s3 = x3(:,l)+j*x3(:,2); 
s4 = x4(:,l)+j*x4(:,2); 

% clear xi's 
clear x 1 ; 
clear x2; 
clear x3; 
clear x4; 

% close all files 
fclose('aH'); 
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<J^******** *********************************** *********************** 
<^****************************************************************** 
% 

% Name: David A Streight 

% 

% Naval Postgraduate School - Monterey, California 

% 

% ssca.m 
% 

% Operating Environment 
% Windows 95 

% 

% Description 

% Calculates the auto spectral density function using the strip 
% spectral correlation algorithm (ssca) 

% 

% Date of last revision 
% 18 November 1996 

% 

% Inputs 
% x - SOI 

% 

% Outputs 

% SDF - the results of the ssca computation. Must be converted 
% and plotted by plotsdf.m 

% 



x = [sl(l:2 A 12,l)].'; 

% prepare variables 
s = size(x); 

N = s( 1 ,2); 

Np = 8; 

xp = [x zeros(l,Np)]; 
w = (hamming(Np))'; 
W = (hamming(N))'; 

L = Np/4; 

P = floor(N/L); 
p = zeros( 1 ,floor(P/L)); 
xT = zeros(l,Np); 
k = -Np/2:Np/2-l; 
dm = zeros(l,Np); 

XTi = zeros(N,Np); 
SXT = zeros(Np,N); 



% number of points in vector 
% course FFT number 

% zero pad input variable for Np length slices 
% hamming window for course FFT 
% hamming window for P point FFT at the end 
% decimation factor 
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% find complex demodulates 
for p = 0:P-1 

xT = fft((xp(l,p*L+l:p*L+Np) .* w), Np); % take Np point FFT of windowed data 

xT = [xT(l,5:8) xT( 1,1:4)]; % fold fft to account for neg freq's 

dm = exp(-j*2*pi*k*p*L/Np); % calculate freq shifts for this set 

XTi(p*L+l,:) = xT .* dm; % freq shift and populate odd col’s of interp 

matrix 
end 



% interpolate ***works for Np = 8 only at this point*** 
for v = 1:P-1 

XTi(L*v,:) = (XTi(L*v-l,:) + XTi(L*v+l,:))./2; % fill in even col's linearly 

end 

XTi(N,:) = XTi(l,:); % hold next to last col for last col values 



%finish up 

xm = [x;x;x;x;x;x;x;x]'; 

%Wm = [W ;W; W; W ;W ; W ; W ; W] ; 
window 

SXT = fft(XTi .* xm, N); 

SDF = SXT.'; 

SDF = [SDF(:, 1 :P-1) SDF(:,P:N)]; 



% matrix of original signal 
% matrix of column-wise hamming 

% take N-point FFT of products 
% transpose after FIT 
% fold fft to account for negative freq's 
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*4* 4* *4^ •A' 4* «df 4* *>i* 4* 'i* 4* ^ *£■* ^ 4* 4; *Jc* !iif 'i' 4* 4* 4* >l-» 4* •Jrf 4* , A‘4 , 'l'»i <,1 i'4 , 'l‘4' ^ >i/ vU vt* nL* vl-» vU vL* vU *1* «I# »i> 4* si/ *1, vL* vL. 

< T* “T 'T' # T' ^ ‘T* -T* *T* -T % - T k ^ - T % ^ ^ ^ 'V ‘T % 'T' ‘T k *T* 'r ^ 'r T* 'T 'T 'T 'T ^ ^ *T* -T* *T“ 'T' -T* ^'V‘T' < T' ,l T' ; r'T''T* ^ ^ ^ ^ Jf* /fv ^ /Js ^ ^ 

{J^****************************************************************** 

% 

% Name: David A Streight 

% 

% Naval Postgraduate School - Monterey, California 

% 

% PLOTSDF: Convert and plot SSCA output 

% 

% Operating Environment 
% Windows 95 

% 

% Description 

% 

% 

% 

% 

% Date of last revision r 

• % 15 November 1996 

% 

% Inputs 

% Matrix of Spectral Density Function values from the SSCA algorithm 
% SDF 

% 

% All variables used in the SSCA computation (no clear executed) 

% 

% Outputs 

% Plot in 3D of the SDF function drawn along lines of cycle frequency 

% 



% determine the span of alpha 
alphamax = 2*N - 2*N/Np - 2; 
alphamin = -2*N; 
alphaspan = alphamax - alphamin; 

% setup intermediate matrices 
I = zeros(alphaspan+l,Np); 

X = zeros(alphaspan+l,Np); 

Y = zeros(alphaspan+l,Np); 

% 

maxalpha = zeros(l,alphamax+l); 

% populate matrices from SDF 
for q = -N/2:N/2-l 
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for k = -Np/2:Np/2 - 1 

I(2*N*k/Np+2*q+2*N+ 1 ,k+Np/2+ 1 ) = abs(SDF(k+Np/2+l,q+N/2+l)); 
end 

end 

% normalize magnitudes 
I = I/max(max(I)); 

% populate f and alpha coordinate matrices 

j = 1:8; 

y = ones(l,8); 

for i = l:alphaspan+l 

X(i,:) = 2*N*(j-Np/2-l)/Np-(i-2*N-l)/2; 

Y(i,:) = ((i-2*N-l).*y)/2; 



end 

% plot the positive alpha half of the plane 
figure(l); 

surf(X(alphaspan-alphamax:alphaspan+l,:),Y(alphaspan- 
alphamax:alphaspan+ 1 ,:),I(alphaspan-alphamax:alphaspan+l 
view(-105,60),grid,axis([-20000 10000 0 8000 0 1]) 

% now plot alpha vs mag of SDF for those values 
for 1 = l:alphamax+l 

maxa(l) = max(I(alphaspan-alphamax+l-l,:)); 

end 

maxalphas = nonzeros(maxa); 
maxalphas = maxalphas'; 
sz = size(maxalphas); 
alphanum = 0:sz( 1 ,2)- 1 ; 

•figure(2); 

plot(alphanum, maxalphas) 

clear i; 
clear j; 
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i^************;}:*;}:* ********************************* ***************** 
%****************************************************************** 



% Naval Postgraduate School - Monterey, California 

% 

% SPECCOA - Cyclostationary Time Difference of Arrival Algorithm 

% 

% Operating Environment 
% Windows 95 

% 

% Description 

% Produces an estimate of the time difference of arrival of a 
% signal of interest between two spatially separated receivers 
% in terms of multiples of the sampling frequency 

% 

% Date of last revision 
% 06 January 1997 

% 

% Inputs 

% p - up sample integer 
% q - down sample integer 

% terms - number of terms used in linear interpolation 
% N - number of samples of decimated sequence to use 
% a - cycle frequency of interest in terms of decimated fs/N 
% SOU - signal captured at receiver 1 (N samples) 

% SOI2 - signal captured at receiver 2 (N samples) 

% *both SOI(s) sampled at rate fs samples per sec 
% 

% Outputs 

% Plot of SPECCOA function versus samples. Peak is TDOA estimate 



sfc 2$c j|c j|c j|c jfc jjc jfc jJc sfc sfc ^ ^ jfc sjc jjc jfc sjc jfc )Jc jjc jjc ^ * jJc jfj jfc ijc jJc jfc ijc )jc ^ jjc ijc )jc )jc * j|c )Jc )jc jfi jJc >fc sfc *Jc jJc 



% 



% Name: David A Streight 



% 



% 



loaddata; 



% load decimated signals 



% resample variables 



p = 256; 

q = i; 



% up sample 
% down sample 

% terms in the linear interpolation 
% total interpolation 



terms = 3; 

I = p/q; 



% cyclostationary variables 
a = 2049 * I; 



% cycle frequency at new sample rate 
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N = 4096 * I; 
fs = 9765.625; 


% length of sequence at new sample rate 


% prepare loop variables 
SXT = zeros(l,N); 
SYXT = zeros(l,N); 


% SCF of SOI1 at alphaO 
% Cross SCF at alphaO 


Sigl =sl(l:(N/I),l).’; 
Sig2 = s2(l:(N/I),l).'; 


% set up vectors 


5011 = resample(Sigl,p,q, terms); 

5012 = resample(Sig2,p,q, terms); 


% resample and set as SOU 
% resample and set as SOI2 


Xf = fft(SOIl,N); 
XT = fftshift(Xf); 
Yf = fft(SOI2,N); 
YT = fftshift(Yf); 


% freq spectrum of SOU 
% adjust for MATLAB fft 
% freq spectrum of SOI2 
% adjust as above 



XT = [zeros( 1 ,a/2) XT zeros(l,a/2)]; % zero pad at each end for correlation loop 
YT = [zeros( 1 ,a/2) YT zeros(l,a/2)]; % ditto 

SXT(1,:) = XT(1,1:N) * conj(XT(l,a+l:N+a)); 

SYXT(1,:) = YT(1,1:N) * conj(XT(l,a+l:N+a)); 

isp = ifft((conj(SXT) .* SYXT), N); % IFFT of complex product from above 



s = -N/2:N/2-l; 
kernel = exp(j*pi*a*s/N); 
sp = fftshift(isp) .* kernel; 
plot(s,real(sp)) 


% time in samples 
% freq kernel at alphaO 
% adjust and multiply 
% plot to find TDOA 
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^^H?**************************************************************** 

% 

% Name: David A Streight 

% 

% Naval Postgraduate School - Monterey, California 

% 

% Closed Form TDOA Geolocation Solution 

% 

% Operating Environment 
% Windows 95 

% 

% Description 

% Produces two WGS84 based 3-dimensional geolocation solutions for a transmitter 
% given the location of each of four observers and the TDOAs for their observations 
% Data outside the algorithm must be used to eliminate the ambiguity. 

% 

% Date of last revision 
% 10 October 1996 

% 

% Inputs 

% WGS84 based x, y and z coordinates for each of four observers 
% xi, yi zi, i = 1,2, 3,4 

% 

% TDOAs for the observers 

% Dij, i = 1 ,2,3,4 and j = 1 ,2,3,4 where i 1 j and Dij = -Dji 

% 

% Outputs 

% WGS84 based x, y and z coordinates for the transmitter 
% rt 1 and rt2 (vectors) 

% 



% positions 

xt = -741941 .47; yt = -5462120.41; zt = 3198242.728; 
xl = -741203.52; yl = -5456323.70; zl = 3208396.55; 
x2 = -735740.31; y2 = -5467456.78; z2 = 3190514.57; 
x3 = -754177.92; y3 = -5459066.31; z3 = 3200764.99; 
x4 = -731248.92; y4 = -5463237.87; z4 = 3198800.64; 



% set up inputs and utility vectors 

tl =0; 

baseline 

t2 = -D12; 

t3 = -D13; 



% relative times of arrival using observer 1 as 

% D12 is defined as tl - t2 thus t2 = -D12 if tl =0 
% ditto 
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t4 = -D14; 



% ditto 



rl = [xl;yl;zl]; 



% observer position vectors 



r2 = [x2;y2;z2]; 
r3 = [x3;y3;z3]; 
r4 = [x4;y4;z4]; 

Euclid_sq 1 = x 1 A 2 + y 1 A 2 + z 1 A 2; % squares of Euclidean distances 

Euclid_sq2 = x2 A 2 + y2 A 2 + z2 A 2; 

Euclid_sq3 = x3 A 2 + y3 A 2 + z3 A 2; 

Euclid_sq4 = x4 A 2 + y4 A 2 + z4 A 2; 

% set up solution matrices 

A = [xl-x2 yl-y2 zl-z2;x2-x3 y2-y3 z2-z3;x3-x4 y3-y4 z3-z4]; 
q = [tl-t2;t2-t3;t3-t4]; 
s = zeros(3,l); 

s(l,l) = (Euclid_sql - Euclid__sq2 - tl A 2 + t2 A 2).*0.5; 
s(2,l) = (Euclid_sq2 - Euclid_sq3 - t2 A 2 + t3 A 2).*0.5; 
s(3, 1) = (Euclid_sq3 - Euclid_sq4 - t3 A 2 + t4 A 2).*0.5; 

g = A\q; % use matrix division (MATLAB)... 



h = A\s; 



% in lieu of the inverse function... 
% to solve this part 



% find terms for at A 2 + bt + c 



a = ((g(l,l)) A 2 + (g(2,l)) A 2 + (g(3,l)) A 2) - 1;% find quadratic term 



b = (g.' * h) + (g.' * r 1 ) + 1 1 ; 



% find first order term 



c = (norm(h - rl). A 2) - tl A 2; 



% find constant 



p = [a 2*b c]; 



% form polynomial in t 



% find the two times as roots of this equation 



t = roots(p); 



% find roots of at A 2 + 2bt + c 



% substitute into t-parametrized equation for r and solve for real part of two solutions 
rtl = real(g .* t(l,l) + h); 
rt2 = real(g .* t(2,l) + h); 
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APPENDIX B. PLOTS FOR ARL:UT TEST DATA 



The following 150 figures are the result of testing 25 the cyclostationary TDOA 
processor for 25 time periods. For each time period, 6 figures appear: one SCF plot for 
each of the four observers grouped two to a figure, one magnitude of SCF vs cycle 
frequency for each of the four observers grouped two to a figure and three SPECCOA 
plots grouped two and one to a figure. 

For the SCF plots, only the positive cycle frequency values are plotted. The 
spectral frequency axis runs from the upper left to the lower right; the cycle frequency 
axis runs from the upper right to the lower left. Both the spectral frequency and the 
cycle frequency are plotted in terms of f s /N where f s = 9765.625 Hz and N = 4096. The 
magnitude is of course the axis pointing to the top of the page. The magnitude of the 
plots is nomalized by the maximum value of the SCF and thus ranges from 0 to 1. 

The magnitude of the SCF versus cycle frequency plots also display only the 
positive cycle frequencies. The horizontal axis is cycle frequency in terms of f s /N as 
above. The vertical axis is again normalized magnitude. These plots merely represent 
looking at the SCF plots from the cycle frequency axis side. They aided in choosing the 
best cycle frequency for SPECCOA. 

Finally, the SPECCOA plots show the output of the speccoa code. The 
horizontal axis represents the TDOA in terms of the numbers of samples of delay. The 
horizontal axis is once again magnitude. 
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Figure B-l - SCF for master (top) and slave #1 (bottom), time 409400 
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Figure B-2 - SCF for slave #2 (top) and slave #3 (bottom), time 409400 
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Figure B-3 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409400 
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Figure B-4 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409400 
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Figure B-5 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409400 
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Figure B-6 - SPECCOA for master and slave #3, time 409400 
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Figure B-7 - SCF for master (top) and slave #1 (bottom), time 409405 
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Figure B-8 - SCF for slave #2 (top) and slave #3 (bottom), time 409405 
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Figure B-9 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409405 
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Figure B-10 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409405 
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Figure B-12 - SPECCOA for master and slave #3, time 409405 
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Figure B-13 - SCF for master (top) and slave #1 (bottom), time 409410 
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Figure B-14 - SCF for slave #2 (top) and slave #3 (bottom), time 409410 
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Figure B-15 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409410 
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Figure B-16 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409410 
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Figure B- 1 7 - SPECCOA for master and slave # 1 (top) and master and slave #2 (bottom), time 4094 1 0 
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Figure B-18 - SPECCOA for master and slave #3, time 409410 
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Figure B-19 - SCF for master (top) and slave #1 (bottom), time 409415 
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Figure B-20 - SCF for slave #2 (top) and slave #3 (bottom), time 409415 
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Figure B-21 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409415 
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Figure B-22 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409415 
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Figure B-23 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409415 
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Figure B-24 - SPECCOA for master and slave #3, time 409415 
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Figure B-25 - SCF for master (top) and slave #1 (bottom), time 409420 
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Figure B-26 - SCF for slave #2 (top) and slave #3 (bottom), time 409420 
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Figure B-27 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409420 
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Figure B-28 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409420 
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Figure B-30 - SPECCOA for master and slave #3, time 409420 
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Figure B-31 - SCF for master (top) and slave #1 (bottom), time 409425 
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Figure B-32 - SCF for slave #2 (top) and slave #3 (bottom), time 409425 
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Figure B-33 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409425 
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Figure B-34 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409425 
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Figure B-35 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409425 
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Figure B-36 - SPECCOA for master and slave #3, time 409425 
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Figure B-37 - SCF for master (top) and slave #1 (bottom), time 409430 
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Figure B-38 - SCF for slave #2 (top) and slave #3 (bottom), time 409430 
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Figure B-39 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409430 
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Figure B-40 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409430 
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Figure B-4 1 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409430 
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Figure B-42 - SPECCOA for master and slave #3, time 409430 
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Figure B-43 - SCF for master (top) and slave #1 (bottom), time 409435 
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Figure B-44 - SCF for slave #2 (top) and slave #3 (bottom), time 409435 
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Figure B-45 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409435 
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Figure B-46 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409435 
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Figure B-50 - SCF for slave #2 (top) and slave #3 (bottom), time 409455 
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Figure B-51 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409455 
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Figure B-52 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409455 
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Figure B-54 - SPECCOA for master and slave #3, time 409455 
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Figure B-55 - SCF for master (top) and slave #1 (bottom), time 409460 
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Figure B-56 - SCF for slave #2 (top) and slave #3 (bottom), time 409460 



0 

5000 



169 





Figure B-57 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409460 
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Figure B-58 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409460 
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Figure B-59 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409460 
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Figure B-60 - SPECCOA for master and slave #3, time 409460 
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Figure B-61 - SCF for master (top) and slave #1 (bottom), time 409465 
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Figure B-62 - SCF for slave #2 (top) and slave #3 (bottom), time 409465 
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Figure B-63 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409465 
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Figure B-64 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409465 
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Figure B-65 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409465 
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Figure B-66 - SPECCOA for master and slave #3, time 409465 
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Figure B-67 - SCF for master (top) and slave #1 (bottom), time 409470 
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Figure B-68 - SCF for slave #2 (top) and slave #3 (bottom), time 409470 
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Figure B-69 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409470 
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Figure B-70 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409470 
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Figure B-72 - SPECCOA for master and slave #3, time 409470 
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Figure B-73 - SCF for master (top) and slave #1 (bottom), time 409475 
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Figure B-74 - SCF for slave #2 (top) and slave #3 (bottom), time 409475 
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Figure B-75 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409475 
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Figure B-76 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409475 
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Figure B-77 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409475 
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Figure B-78 - SPECCOA for master and slave #3, time 409475 
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Figure B-79 - SCF for master (top) and slave #1 (bottom), time 409480 
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Figure B-80 - SCF for slave #2 (top) and slave #3 (bottom), time 409480 
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Figure B-81 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409480 
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Figure B-82 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409480 
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Figure B-83 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409480 
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Figure B-84 - SPECCOA for master and slave #3, time 409480 
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-85 - SCF for master (top) and slave #1 (bottom). 
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Figure B-86 - SCF for slave #2 (top) and slave #3 (bottom), time 409485 
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Figure B-87 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409485 
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Figure B-88 - Cycle frcq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409485 
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Figure B-89 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409485 
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Figure B-90 - SPECCOA for master and slave #3, time 409485 
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Figure B-91 - SCF for master (top) and slave #1 (bottom), time 409490 
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Figure B-92 - SCF for slave #2 (top) and slave #3 (bottom), time 409490 
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Figure B-93 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409490 
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Figure B-94 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409490 
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Figure B-95 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409490 
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Figure B-96 - SPECCOA for master and slave #3, time 409490 



209 




1000 



2000 



4000 



3000 



0 

5000 




1000 



2000 



4000 



3000 



5000 



Figure B-97 - SCF for master (top) and slave #1 (bottom), time 409495 
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Figure B-98 - SCF for slave #2 (top) and slave #3 (bottom), time 409495 
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Figure B-99 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409495 
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Figure B-100 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409495 
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Figure B-101 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409495 
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Figure B-102 - SPECCOA for master and slave #3, time 409495 
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Figure B-103 - SCF for master (top) and slave #1 (bottom), time 409500 
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Figure B-104 - SCF for slave #2 (top) and slave #3 (bottom), time 409500 
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Figure B-105 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409500 
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Figure B-106 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409500 
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Figure B-107 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409505 
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Figure B-108 - SPECCOA for master and slave #3, time 409505 
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Figure B-109 - SCF for master (top) and slave #1 (bottom), time 409505 
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Figure B- 1 1 1 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409505 
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Figure B-l 12 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409505 
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Figure B- 113- SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409505 
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Figure B-l 14 - SPECCOA for master and slave #3, time 409505 
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Figure B-l 15 - SCF for master (top) and slave #1 (bottom), time 409510 
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Figure B- 116- SCF for slave #2 (top) and slave #3 (bottom), time 409510 
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Figure B- 117 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409510 
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Figure B- 118- Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409510 
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Figure B-l 19 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409510 
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Figure B-120 - SPECCOA for master and slave #3, time 409510 



233 




1000 



2000 



4000 



3000 



5000 




3000 



1000 



2000 



0 

5000 



-5000 



4000 



Figure B- 121 - SCF for master (top) and slave #1 (bottom), time 409515 
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Figure B-122 - SCF for slave #2 (top) and slave #3 (bottom), time 409515 
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Figure B-123 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409515 
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Figure B-124 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409515 
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Figure B-126 - SPECCOA for master and slave #3, time 409515 
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Figure B-127 - SCF for master (top) and slave #1 (bottom), time 409520 



240 




Figure B-128 - SCF for slave #2 (top) and slave #3 (bottom), time 409520 
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Figure B-129 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409520 
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Figure B-I30 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409520 
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Figure B- 131 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409520 
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Figure B-132 - SPECCOA for master and slave #3, time 409520 
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Figure B-133 - SCF for master (top) and slave #1 (bottom), time 409525 
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Figure B-134 
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Figure B-135 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409525 
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Figure B-136 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409525 
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Figure B-137 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409525 
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Figure B-138 - SPECCOA for master and slave #3, time 409525 
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Figure B-139 - SCF for master (top) and slave #1 (bottom), time 409530 
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Figure B-140 - SCF for slave #2 (top) and slave #3 (bottom), time 409530 
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Figure B-141 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409530 
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Figure B-142 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409530 
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Figure B-143 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409530 
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Figure B-144 - SPECCOA for master and slave #3, time 409530 
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Figure B-145 - SCF for master (top) and slave #1 (bottom), time 409535 
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Figure B-146 - 
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Figure B-147 - Cycle freq vs Max-magnitude of SCF for master (top) and slave #1 (bottom), time 409535 
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Figure B-148 - Cycle freq vs Max-magnitude of SCF for slave #2 (top) and slave #3 (bottom), time 409535 
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Figure B-149 - SPECCOA for master and slave #1 (top) and master and slave #2 (bottom), time 409535 
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Figure B-150 - SPECCOA for master and slave #3, time 409535 
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